Numerical Analysis

PART1  <求解方程>

1，二分法

def bisect(f,a,b,TOL=0.000004):
    u_a = a
    u_b = b
    while(u_b-u_a)/2.0 > TOL:
        c = (u_a+u_b)/2.0
        if f(c) == 0:
            break
        if f(u_a)*f(c) < 0:
            u_b = c
        else:
            u_a = c

    u_c = (u_a + u_b) / 2.0
    return u_c

f = lambda x: x*x*x + x - 1
ret = bisect(f,-1.0,1.0)
print(ret)

print(f(ret))

        
    

2,不动点迭代法求解方程(FPI)

例如求cosx - x = 0 方程 事实上可以化成x = cosx

x^3 + x - 1 = 0 也可以化成 x=1 - x^3

等式左边是x 右边设为f(x)

即是:x = f(x)

过程如下:
x1 = f(x0)
x2 = f(x1)
x3 = f(x2)
........

例如求cosx = x  即(cosx -  x = 0）的解:  他的不动点是0.7390851332

def fpi(f,x0,k):
    xvalues = []
    xvalues.append(x0)
    for i in range(k):
        xvalues.append(f(xvalues[i]))
    print(xvalues)
    return xvalues[-1]  # return last

import math
f = lambda x:math.cos(x)
v = fpi(f,0,10)
print(v)

上述代码迭代十次 可以看到生成的数值为:

[0, 1.0, 0.5403023058681398, 0.8575532158463933, 0.6542897904977792, 0.7934803587425655, 0.7013687736227566, 
0.7639596829006542, 0.7221024250267077, 0.7504177617637605, 0.7314040424225098, 0.7442373549005569, 
0.7356047404363473, 0.7414250866101093, 0.7375068905132428, 0.7401473355678757, 0.7383692041223232, 
0.739567202212256, 0.7387603198742114, 0.7393038923969057, 0.7389377567153446]

20次基本达到了与精确解

PART2  <求解方程组>

3，雅克比(jacobi)和高斯-赛德尔(Gauss-Seidel) 解方程组

结论:高斯-赛德尔的收敛速度比雅克比快的多。Gauss-Seidel20次迭代达到精确解，jacobi 20次还没得到

还有个高斯-赛德尔SOR 实现起来比较简单。

D代表对角阵，除了对角元素其他都为0

L,除了对角元素往下的元素，其他都为0

U,除了对角元素往上的元素，其他都为0

注意这里的LU和 LU分解是有本质的区别。

高斯赛的尔 用的我蓝色标识的字:

import numpy as np

class GMatrix:
    @staticmethod
    def Invert(matrix):
        nm = np.linalg.inv(matrix.mmatrix)
        gm = GMatrix()
        gm.assignNumPyMatrix(nm)
        return gm



    def __init__(self, nx=0, ny=0):
        self.rows = nx
        self.columns = ny
        self.mmatrix = np.matrix(np.zeros((nx, ny)))

    def assignNumPyMatrix(self, matrix):
        self.mmatrix = matrix

    def changeShape(self, nx, ny):
        self.mmatrix = np.matrix(np.zeros((nx, ny)))

    def debug(self):
        print("{
", ">>SHAPE:", (self.mmatrix.shape))
        print(self.mmatrix)
        print('}
')

    def numPyMatrix(self):
        return self.mmatrix

    def set(self, i, j, var):
        self.mmatrix[i - 1, j - 1] = var

    def get(self, i, j):
        return self.mmatrix[i - 1, j - 1]

    def L(self):
        n = self.rows
        m = GMatrix(n, n)

        for j in range(1, n, 1):
            for i in range(j + 1, n + 1):
                # print(i, j)
                m.set(i, j, self.get(i, j))
        return m

    def U(self):
        n = self.rows
        m = GMatrix(n, n)
        for i in range(1, n + 1, 1):
            for j in range(i + 1, n + 1):
                m.set(i, j, self.get(i, j))
        return m

    def D(self):
        """
        this matrix diagonal values
        1,make a empty matrix
        2,copy diagonal value to this empty matrix
        """
        n = self.rows
        m = GMatrix(n, n)
        for i in range(self.rows):
            index = i + 1
            m.set(index, index, self.get(index, index))
        return m

    def DInvert(self):
        """
        invert diagonal matrix !
        :return:
        """
        n = self.rows
        m = GMatrix(n, n)
        for i in range(n):
            index = i + 1
            m.set(index, index, 1 / self.get(index, index))
        return m

    def __mul__(self, other):
        nm = self.mmatrix * other.mmatrix  # USE NumPy matrix * matrix
        m = GMatrix()
        m.assignNumPyMatrix(nm)
        return m

    def __add__(self, other):
        nm = self.mmatrix + other.mmatrix  # USE NumPy matrix + matrix
        m = GMatrix()
        m.assignNumPyMatrix(nm)
        return m

    def __sub__(self, other):
        nm = self.mmatrix - other.mmatrix  # USE NumPy matrix - matrix
        m = GMatrix()
        m.assignNumPyMatrix(nm)
        return m


def Jacobi_Solver(x0, L, U, invertD, b,iter =20):
    """
    x0   = INIT_VECTOR_VALUE
    xk+1 = Inv(D)*[ b-(L+U)xk ]  ,k = 0,1,2,......
    :return:
    """
    nextX = x0
    for x in range(iter):
        nextX = invertD * (b - (L + U) * nextX)
    return nextX


def Gauss_Seidel_Solver(x0, L, U, D, b,iter =10):
    """
    xk+1 = inv(L+D) * (b- U*xk)

    """
    nextX = x0
    for x in range(iter):
        nextX = GMatrix.Invert(L+D) * (b - U * nextX)
    return nextX




def unit_test_matrix():
    m = GMatrix(2, 2)
    m.set(1, 1, 3)
    m.set(1, 2, 1)
    m.set(2, 1, 1)
    m.set(2, 2, 2)

    # and also can this
    # m.numPyMatrix()[0, 0] = 3
    # m.numPyMatrix()[0, 1] = 1
    # m.numPyMatrix()[1, 0] = 1
    # m.numPyMatrix()[1, 1] = 2
    m.debug()
    m.D().debug()
    m.L().debug()
    m.U().debug()
    m.DInvert().debug()


def unit_test_jacobi():
    """
    |  3 1  |   |u|      |5|
    |       |   | |   =  | |
    |  1 2  |   |v|      |5|

    --------
    AX = b
    --------
    :return:
    """

    # CONSTRUCT THE A MATRIX
    A = GMatrix(2, 2)
    A.set(1, 1, 3)  # A11 = 3
    A.set(1, 2, 1)  # A12 = 1
    A.set(2, 1, 1)  # A21 = 1
    A.set(2, 2, 2)  # A22 = 2
    #A.debug()

    # CONSTRUCT THE b vector
    b = GMatrix(2, 1)
    b.set(1, 1, 5)
    b.set(2, 1, 5)

    # INIT zero vector to iter
    x0 = GMatrix(2,1) # [0,0]T vector initialize

    invert_D = A.DInvert()

    L = A.L()
    U = A.U()

    print("get the jacobi result ")
    Jacobi_Solver(x0,L,U,invert_D,b).debug()


def unit_test_Gauss_Seidel():
    """
        |  3 1  |   |u|      |5|
        |       |   | |   =  | |
        |  1 2  |   |v|      |5|

        --------
        AX = b
        --------
        :return:
        """

    # CONSTRUCT THE A MATRIX
    A = GMatrix(2, 2)
    A.set(1, 1, 3)  # A11 = 3
    A.set(1, 2, 1)  # A12 = 1
    A.set(2, 1, 1)  # A21 = 1
    A.set(2, 2, 2)  # A22 = 2


    # CONSTRUCT THE b vector
    b = GMatrix(2, 1)
    b.set(1, 1, 5)
    b.set(2, 1, 5)

    # INIT zero vector to iter
    x0 = GMatrix(2, 1)  # [0,0]T vector initialize

    D = A.D()
    L = A.L()
    U = A.U()

    print("get the Gauss_Seidel result ")
    Gauss_Seidel_Solver(x0, L, U, D, b).debug()

unit_test_jacobi()
unit_test_Gauss_Seidel()

4,楚列斯基分解choleskyDecomposition

这类情况只适合于正定矩阵。

楚列斯基分解的目的就是将正定矩阵分解成A = transpose(R) * R .

求解方程组直接使用回代法即可求出向量解,

这类方法不是迭代法，最重要的是分解矩阵.

所以是要求出来这个上三角矩阵R,比如分解如下矩阵：

[[  4.  -2.   2.]
 [ -2.   2.  -4.]
 [  2.  -4.  11.]]

import numpy as np
import math


class GMatrix:
    def __init__(self, shape):
        self.nn = shape
        self.Mat = np.zeros((shape, shape))

    def setRawData(self, *args):
        if len(args) != self.nn * self.nn:
            raise ValueError

        deserialize = []
        for i in range(self.nn):
            rhtmp = []
            for j in range(self.nn):
                rhtmp.append(args[i * self.nn + j])
            deserialize.append(rhtmp)
        self.tupleToMatrix(deserialize)

    def tupleToMatrix(self, tupleObj):
        """
        :param tupleObj: ([1,2],[3,4])
        will set matrix above tuple obj
        """
        for i in range(self.nn):
            # get row data of args
            row_data = tupleObj[i]
            for j in range(self.nn):
                self.Mat[i, j] = row_data[j]

    def setRowsData(self, *args):
        if len(args) != self.nn:
            raise ValueError
        self.tupleToMatrix(args)

    def __str__(self):
        return str(self.Mat)

    def divideConstant(self, var):
        self.Mat[:, :] = self.Mat[:, :] / var

    def __add__(self, other):
        self.Mat[:, :] = self.Mat[:, :] + other.mat[:, :]

    def __sub__(self, other):
        self.Mat[:, :] = self.Mat[:, :] - other.mat[:, :]

    def nextLayer(self):
        array = self.Mat[1:, 1:]
        m = GMatrix(a.Mat.shape[0])
        m.Mat = array
        return m

    def getupu(self):
        self.Mat = np.triu(self.Mat)
        return self

def choleskyDecomposition(op=GMatrix(2)):
    """ DEFINE MATRIX:
    R11 R12 ... R1n
    R21 R22 ... R2n
     .   .       .
     .   .       .
    Rn1 Rn2 ... Rnn

    """
    if len(op.Mat.flatten()) == 1:
        op.Mat[op.nn - 1:, op.nn - 1:] = np.sqrt(op.Mat[op.nn - 1:, op.nn - 1:])
        return

    # First make the matrix A11 to sqrt(A11)
    op.Mat[0, 0] = math.sqrt(op.Mat[0, 0])  # R11
    op.Mat[0, 1:] = op.Mat[0, 1:] / op.Mat[0,0]  # (R12 R13 .... R1n) = (R12 R13 ...R1n) / R11

    # cal u.transpose(u)
    u   = np.matrix(op.Mat[0,1:])      # it's a row vector
    ut  = np.transpose(u)              # it's a column vector

    uut = ut*u
    # extract_layer - uut
    op.Mat[1:,1:] = op.Mat[1:,1:] - uut
    choleskyDecomposition(op.nextLayer())



a = GMatrix(3)

a.setRawData(4, -2, 2, -2, 2, -4, 2, -4, 11)
print(a)
choleskyDecomposition(a)
print(a.getupu())

分解结果R:

[[ 2. -1.  1.]
 [ 0.  1. -3.]
 [ 0.  0.  1.]]

Mathematica 检查结果:

 

5,共轭梯度法:

6,预条件共轭梯度法:

BVP:

差分法比较简单略.

排列法求解BVP问题:

如下例，n=2书中已经给答案，n=4的矩阵未给出，顺便自己解下矩阵n=4的情况。

算出c1,c2,c3,c4:

PDE 解热方程:

这次将时间方向放入坐标轴，可以用三维图像直接看到 在时间步上，热方程是如何传递的.

 

clc;
clear w a b c d M N f l r m sigma ;
xl = 0;
xr = 1;
yb = 0;
yt = 1;

M = 10;
N = 250;

f = @(x) sin(2*pi*x).^2;
l = @(t) 0*t;
r = @(t) 0*t;

D =1;
h = (xr-xl) / M
k = (yt-yb) / N
m = M -1;
n=N;
sigma = D*k/(h*h);

lside = l(yb+(0:n)*k);
rside = r(yb+(0:n)*k);



% 定义矩阵a
a = diag(1-2*sigma*ones(m,1));
a = a + diag(sigma*ones(m-1,1),1);  %自身加上  往右offset 1列 sigma倍数
a = a + diag(sigma*ones(m-1,1),-1); %自身加上  往左偏移 1列sigma倍数
a

% 设置w第一列 ,初值条件,注意 要转置，因为 设置 w的第一列向量
w(:,1) = f(xl + (1:m) * h)'

disp('START FOR LOOP')
% 给w 在时间250 迭代成列向量，则w是251列向量
for j = 1:n
    rhs = w(:,1)';
    w(:,j+1) = a*w(:,j) ;
end
disp('END FOR LOOP')

w = [lside;w;rside];
x = (0:m+1)*h; t= (0:n)*k;
mesh(x,t,w')
view(60,30);
axis([xl xr yb yt -1 1]);