# DickLiTQ / NumAnalysis

This repository is used to create a project on lesson "Numerical Analysis"
Python
Fetching latest commit…
Cannot retrieve the latest commit at this time.
Type Name Latest commit message Commit time
Failed to load latest commit information.
Figure_Interpolation&Fitting.png
Figure_ODE.png
LU.py
LUDecomposition.gif
ODE.gif
ODE_eg1.gif
__init__.py
differential.py
example_of_numerical_integral.png
integral.py
interpolation.py
linrecursive.gif
nonlineareq.py

## 数值分析/NumAnalysis

This repository is used to create a project on lesson "Numerical Analysis" and practically use the method we have learned in problem solving. As a learner who knows only a little about computer and programming, I know that there must exist difficulties in the transformation from theory to practice. It takes long time when the theory generally improves efficiency as technical tools.

From our perspective, some methods in this lesson are really interesting and fantastic. Maybe some much advanced technique has been used in popular modules, however, it is not important while we hope that the method in our textbook can solve the problem we faced and adapted to our demand.

• 非线性方程数值解/Numerical Solution to Nonlinear Equation
• 线性方程组LU分解/System of Linear Equation——LU Decomposition
• 线性方程组迭代法/Recursive way to solve the System of Linear Equation
• 插值与拟合/Interpolation and Fitting
• 数值积分/Numerical Integral
• 常微分方程数值解/Numerical Solution in Ordinary Differential Equation

## 非线性方程数值解/Numerical Solution to Nonlinear Equation

### 迭代法

#### Recursive Method

```import numpy as np
import sympy as sp
g = lambda x: x**2-3*x**(1/2)+4 # Here input your recursive function
x0 = 1 # Your initial point
n = 10 # Your iteration times
residual = 0.001 # Your tolerance

def Recursive_Method(g,x0,n,residual):
x = np.zeros(n+1)
x[0] = x0
for index in range(1,n+1):
x[index] = g(x[index-1])
difference = x[index] - x[index-1]
if np.abs(x[index]-x[index-1])<residual:
print("Iteration %r: x=%r, difference=%r <= Residual=%r" %(index,x[index],difference,residual))
break
else:
print("Iteration %r: x=%r, difference=%r" %(index,x[index],difference))
print("Terminal: The final result is %r" % x[index])
return x[index]

Recursive_Method(g,x0,n,residual)```

```g = lambda x: 0.5*(10-x**3)**0.5
Recursive_Method(g,1.5,10,1e-5)```

``````Iteration 1: x=1.2869537676233751, difference=-0.2130462323766249
Iteration 2: x=1.4025408035395783, difference=0.11558703591620323
Iteration 3: x=1.3454583740232942, difference=-0.057082429516284172
Iteration 4: x=1.3751702528160383, difference=0.029711878792744173
Iteration 5: x=1.3600941927617329, difference=-0.015076060054305396
Iteration 6: x=1.3678469675921328, difference=0.0077527748303998223
Iteration 7: x=1.3638870038840212, difference=-0.0039599637081115802
Iteration 8: x=1.3659167333900399, difference=0.0020297295060187626
Iteration 9: x=1.3648782171936771, difference=-0.0010385161963628597
Iteration 10: x=1.3654100611699569, difference=0.00053184397627981106
Terminal: The final result is 1.3654100611699569
``````

1. 本部分尚未加入是否发散的判断，因此要根据迭代结果进行发散的判断
2. 请自行计算迭代式，在输入过程中注意指数与根号的输入规范

### 牛顿法

#### Newton Method

```f = lambda x: x**2+3*x+1 # Your function f(x)
f1 = lambda x: 2*x+3 # First order deviation of f(x)
x0 = 1 # The initial point
n = 100 # Iteration time
residual = 1e-5 # Tolerance

def Newton_Recursive(f,f1,x0,n,residual):
x = np.zeros(n+1)
x[0] = x0
for index in range(1,n+1):
x[index] = x[index-1]-f(x[index-1])/f1(x[index-1])
difference = x[index] - x[index-1]
if np.abs(x[index]-x[index-1])<residual:
print("Iteration %r: x=%r, difference=%r <= Residual=%r" %(index,x[index],difference,residual))
break
else:
print("Iteration %r: x=%r, difference=%r" %(index,x[index],difference))
print("Terminal: The final result is %r" % x[index])
return x[index]

Newton_Recursive(f,f1,x0,n,residual)```

```f = lambda x: x**3+4*x**2-10
f1 = lambda x: 3*x**2+8*x
Newton_Recursive(f,f1,1.5,10,1e-5)```

``````Iteration 1: x=1.3733333333333333, difference=-0.12666666666666671
Iteration 2: x=1.3652620148746266, difference=-0.0080713184587066777
Iteration 3: x=1.3652300139161466, difference=-3.2000958479994068e-05
Iteration 4: x=1.3652300134140969, difference=-5.0204973511824846e-10 <= Residual=1e-05
Terminal: The final result is 1.3652300134140969
``````

### 牛顿-拉夫逊法

#### Newton-Ralfsnn Method

，其中为Jacobi矩阵（全体一阶偏导矩阵）

```def Jacobi_Matrix(f1,f2,f3,ak,bk,ck):
a = sp.Symbol('a')
b = sp.Symbol('b')
c = sp.Symbol('c')
F = [f1,f2,f3]
X = [a,b,c]
n = len(F)
#    F_1 = np.array([[a,b,c],
#                    [a,b,c],
#                    [a,b,c]])
F_1v = np.zeros((n,n))
for i in range(n):
for j in range(n):
#            F_1[i,j] = sp.diff(F[i],X[j])
storage = sp.diff(F[i],X[j])
F_1v[i,j] = storage.evalf(subs = {a:ak,b:bk,c:ck})
#F1v = F_1.evalf(subs = {X:X0})
return np.linalg.inv(F_1v)

def Newton_Recursive_Multi(f1,f2,f3,a0,b0,c0,n):
a = sp.Symbol('a')
b = sp.Symbol('b')
c = sp.Symbol('c')
x = np.zeros((3,n+1))
x[0,0] = a0
x[1,0] = b0
x[2,0] = c0
for k in range(n):
if np.min(np.abs(x[:,k+1]-x[:,k]))>100:
return "Error with divergence"
#            print("Error")
else:
if np.max(np.abs(x[:,k+1]-x[:,k]))>1e-5:
F1 = f1.evalf(subs = {a:x[0,k],b:x[1,k],c:x[2,k]})
F2 = f2.evalf(subs = {a:x[0,k],b:x[1,k],c:x[2,k]})
F3 = f3.evalf(subs = {a:x[0,k],b:x[1,k],c:x[2,k]})
F1 = float(F1)
F2 = float(F2)
F3 = float(F3)
F = np.array([F1,F2,F3])
x[:,k+1] = x[:,k]-Jacobi_Matrix(f1,f2,f3,x[0,k],x[1,k],x[2,k])@F
else:
return x[:,:k+1]
#                print("Error")
return x```

```# f = [f1,f2,f3,f4,f5]'=0
# a* = 1 b* = 3 c* = 2
a = sp.Symbol('a')
b = sp.Symbol('b')
c = sp.Symbol('c')
f1 = a**2 + b - 2*c
f2 = 3*a + 2**b -c**3 - 3
f3 = a + b + c - 6
a0 = 10
b0 = 25
c0 = 80
n = 100
solution = Newton_Recursive_Multi(f1,f2,f3,a0,b0,c0,n)```

```# f = [f1,f2,f3,f4,f5]'=0
# a* = 1 b* = 3 c* = 2
a = sp.Symbol('a')
b = sp.Symbol('b')
c = sp.Symbol('c')
f1 = a**2 + b - 2*c
f2 = 3*a + 2**b -c**3 - 3
f3 = a + b + c - 6
a0 = 1
b0 = 2.5
c0 = 3.3
n = 100
solution = Newton_Recursive_Multi(f1,f2,f3,a0,b0,c0,n)```

## 线性方程组LU分解/System of Linear Equation——LU Decomposition

### 杜利特尔分解

#### Doolittle Decomposition

```import numpy as np

""" Input Example
A = np.array([[1,3,5,9],
[2,5,2,5],
[9,3,4,1],
[1,10,0,19]]) # Here input your coefficient matrix
b = np.array([10,24,31,42]) # Here input the vector
"""

def check(A):
row = A.shape[0]
col = A.shape[1]
if row!=col:
print("Input error: A is not a square matrix")
return 0
else:
if np.linalg.det(A)==0:
print("The determination of A is equal to zero")
return 0
else:
if row == 1:
print("The dimension of matrix is 1")
return 0
else:
return row```

```def Decomposition(A):
if check(A) == 0:
print("Error")
else:
print("det(A)=%r"%np.linalg.det(A))
dim = check(A)
L = np.eye(dim)
U = np.zeros_like(A)
U[0,:]=A[0,:]
L[1:,0]=A[1:,0]/U[0,0]
for r in range(1,dim):
for l in range(1,r):
L[r,l]=1/U[l,l]*(A[r,l]-L[r,:l]@U[:l,l])
for u in range(r,dim):
U[r,u]=A[r,u]-L[r,:r]@U[:r,u]
print("L=\n",L,"\n","U=\n",U)
return L,U```

```A = np.array([[1,3,5],
[2,5,2],
[9,3,4]
])
Decomposition(A)```

```det(A)=-151.0
L=
[[  1.   0.   0.]
[  2.   1.   0.]
[  9.  24.   1.]]
U=
[[  1   3   5]
[  0  -1  -8]
[  0   0 151]]```

```def Solve(A,b):
L,U = Decomposition(A)
y = np.linalg.solve(L,b)
x = np.linalg.solve(U,y)
print("y=\n",y,"\n","x=\n",x)
return y,x```

## 线性方程组迭代法/Recursive way to solve the System of Linear Equation

### Jacobi迭代

#### Jacobi Recursion

Jacobi迭代的格式为

```import numpy as np
import sympy as sp

def Decomposition(A):
dim = A.shape[0]
D = np.zeros_like(A)
L = np.zeros_like(A)
U = np.zeros_like(A)
for r in range(dim):
D[r,r] = A[r,r]
for c in range(dim):
if c<r:
L[r,c] = -A[r,c]
if c>r:
U[r,c] = -A[r,c]
return D,L,U```

```def Jacobi(A,x0,b,n):
D,L,U = Decomposition(A)
D_1 = np.linalg.inv(D)
print("Inverse of D=\n",D_1)
g = lambda x: D_1@(L+U)@x + D_1@b
x = np.zeros((A.shape[0],n+1))
x[:,0] = x0
for iteration in range(1,n+1):
x[:,iteration] = g(x[:,iteration-1])
print("Iteration %d:"%iteration,"x=",x[:,iteration])
return x```

```A = np.array([[8,-3,2],
[4,11,-1],
[2,1,4]])
b = np.array([20,33,12])
x0 = np.array([0,0,0])
Jacobi(A,x0,b,10)```

``````Inverse of D=
[[ 0.125       0.          0.        ]
[ 0.          0.09090909  0.        ]
[ 0.          0.          0.25      ]]
Iteration 1: x= [ 2.5  3.   3. ]
Iteration 2: x= [ 2.875       2.36363636  1.        ]
Iteration 3: x= [ 3.13636364  2.04545455  0.97159091]
Iteration 4: x= [ 3.02414773  1.94783058  0.92045455]
Iteration 5: x= [ 3.00032283  1.9839876   1.00096849]
Iteration 6: x= [ 2.99375323  1.99997065  1.00384168]
Iteration 7: x= [ 2.99902857  2.0026208   1.00313072]
Iteration 8: x= [ 3.00020012  2.00063786  0.99983051]
Iteration 9: x= [ 3.00028157  1.99991182  0.99974048]
Iteration 10: x= [ 3.00003181  1.99987402  0.99988126]
``````

### Gauss-Seidel迭代

#### Gauss-Seidel Recursion

Gauss-Seidel迭代的格式为 类似地，我们用以下代码进行实现：

```def GS(A,x0,b,n):
D,L,U = Decomposition(A)
D_L = D-L
D_L_1 = np.linalg.inv(D_L)
print("Inverse of D-L=\n",D_L_1)
g = lambda x: D_L_1@U@x + D_L_1@b
x = np.zeros((A.shape[0],n+1))
x[:,0] = x0
for iteration in range(1,n+1):
x[:,iteration] = g(x[:,iteration-1])
print("Iteration %d:"%iteration,"x=",x[:,iteration])
return x```

Gauss-Seidel与Jacobi迭代的区别在于，每次迭代的过程中就代入迭代值加速收敛，因此其收敛速度更快，对于上面方程组，Gauss-Seidel迭代法的10次迭代结果为：

``````Inverse of D-L=
[[ 0.125       0.          0.        ]
[-0.04545455  0.09090909  0.        ]
[-0.05113636 -0.02272727  0.25      ]]
Iteration 1: x= [ 2.5         2.09090909  1.22727273]
Iteration 2: x= [ 2.97727273  2.02892562  1.00413223]
Iteration 3: x= [ 3.00981405  1.99680691  0.99589125]
Iteration 4: x= [ 2.99982978  1.99968838  1.00016302]
Iteration 5: x= [ 2.99984239  2.00007213  1.00006077]
Iteration 6: x= [ 3.00001186  2.00000121  0.99999377]
Iteration 7: x= [ 3.00000201  1.9999987   0.99999932]
Iteration 8: x= [ 2.99999968  2.00000005  1.00000014]
Iteration 9: x= [ 2.99999998  2.00000002  1.        ]
Iteration 10: x= [ 3.00000001  2.          1.        ]
``````

## 插值与拟合/Interpolation and Fitting

### 多项式插值

#### Polynomial Interpolation

```import numpy as np
import sympy as sp
""" Data Input: We input data in an array which has the form of (x,y)

data = np.array([[1,-1],
[2,-1],
[3,1]
])
"""
def Polynomial(data):
dim = data.shape[0]
A = np.zeros((dim,dim))
for c in range(dim):
g = lambda x:x**c
A[:,c] = g(data[:,0])
soln = np.linalg.solve(A,data[:,1])
for index in range(dim):
print("c%d = %r"%(index,soln[index]))
return soln

def Polynomial_predict(data,x):
coefficient = Polynomial(data)
result = 0
for index in range(len(coefficient)):
g = lambda x: x**index
result = result + g(x)*coefficient[index]
return result```

(x1,y1)=(1,0)

(x2,y2)=(2,4)

(x3,y3)=(3,8)

```data = np.array([[1,0],
[2,4],
[3,8]])
Polynomial(data)
Polynomial_predict(data,1.5)```

``````array([-4.,  4., -0.])
2.0
``````

### Lagrange插值

#### Lagrangian Interpolation

```def li_calc(data,i):
dim = data.shape[0]
x = sp.Symbol('x')
li = 1
for index in range(dim):
if index != i:
li = li*(x-data[index,0])/(data[i,0]-data[index,0])
else:
li = li
return li

def Lagrangian(data):
dim = data.shape[0]
L = 0
for index in range(dim):
L = L +li_calc(data,index)*data[index,1]
return L.expand()

def Lagrangian_predict(data,y):
result = Lagrangian(data)
x = sp.Symbol('x')
return result.evalf(subs = {x:y})```

(x1,y1)=(1,-1)

(x2,y2)=(2,-1)

(x3,y3)=(3,1)

```data = np.array([[1,-1],
[2,-1],
[3,1]])
Lagrangian(data)
Lagrangian_predict(data,1.5)```

``````x**2 - 3*x + 1
-1.25000000000000
``````

### Newton插值

#### Newton Interpolation

```def Divided_difference(data):
dim = data.shape[0]
D = np.zeros((dim,dim+1))
D[:,0] = data[:,0]
D[:,1] = data[:,1]
diff = 1
for column in range(2,dim+1):
for row in range(diff,dim):
D[row,column] = (D[row,column-1]-D[row-1,column-1])/(D[row,0]-D[row-diff,0])
diff = diff+1
return D```

```def Newton(data):
Coefficient = Divided_difference(data)
x = sp.Symbol('x')
Newton = 0
X = 1
dim = data.shape[0]
for index in range(dim):
Newton = Newton + Coefficient[index,index+1]*X
X = X*(x-Coefficient[index,0])
return Newton.expand()
def Newton_predict(data,y):
result = Newton(data)
x = sp.Symbol('x')
return result.evalf(subs = {x:y})```

(0,2), (1,-3), (2,-6), (3,11)。

```data = np.array([[0,2],
[1,-3],
[2,-6],
[3,11]])
Divided_difference(data)
Newton(data)
Newton_predict(data,1.5)```

``````array([[  0.,   2.,   0.,   0.,   0.],
[  1.,  -3.,  -5.,   0.,   0.],
[  2.,  -6.,  -3.,   1.,   0.],
[  3.,  11.,  17.,  10.,   3.]])
3.0*x**3 - 8.0*x**2 + 2.0
-5.87500000000000
``````

### 最小二乘拟合

#### Least Square

```def LeastSquare(data,g1,g2,g3):
dim = data.shape[0]
A = np.zeros((dim,3))
A[:,0] = g1(data[:,0])
A[:,1] = g2(data[:,0])
A[:,2] = g3(data[:,0])
if np.linalg.det(A.T@A)!=0 :
result = (np.linalg.inv(A.T@A))@(A.T)@(data[:,1])
print("Coefficient: ",result)
return result
else:
print("The inverse of A'A doesn't exist")
return 0```

```data = np.array([[1,-3],
[2,-6],
[3,11],
[4,7],
[5,14]])
g1 = lambda x: np.exp(x)
g2 = lambda x: x
g3 = lambda x: 1/x
LeastSquare(data,g1,g2,g3)```

``````Coefficient:  [ 0.02822754  2.23012035 -7.19255383]
``````

### 插值与拟合绘图

#### Plot Illustration

```def plot(data):
xmin = np.min(data[:,0])
xmax = np.max(data[:,0])
ymin = np.min(data[:,1])
ymax = np.max(data[:,1])
x0 = np.linspace(xmin-0.5,xmax+0.5,100)
x = sp.Symbol('x')
fig, ax = plt.subplots(dpi=120)
ax.set_xlim(xmin-1,xmax+1)
ax.set_ylim(ymin-1,ymax+1)
ax.scatter(data[:,0],data[:,1],marker='v')
Poly = Polynomial_predict(data,x0)
Lagr = Lagrangian(data)
Newt = Newton(data)
L = np.zeros(100)
N = np.zeros(100)
index = 0
for tick in x0:
L[index] = Lagr.evalf(subs = {x:tick})
N[index] = Newt.evalf(subs = {x:tick})
index = index+1
ax.plot(x0,Poly,'r-',label='Ploynomial',alpha=0.3)
ax.plot(x0,L,'b-',label='Lagrangian',alpha=0.3)
ax.plot(x0,N,'g-',label='Newton',alpha=0.3)
ax.legend()
ax.set_title("Plot")
fig.show()
plot(data)```

## 数值积分/Numerical Integral

### Newton-Cotes公式

#### 梯形公式

##### Trapezium Rule

```def trapezium(data,a,b):
" Integral from ath data to bth data "
#dim = b-a
x = data[a-1,0]
y = data[b-1,0]
fx = data[a-1,1]
fy = data[b-1,1]
return (y-x)/2*(fx+fy)```

#### Sinpson公式

##### Sinpson Rule

```def sinpson(data,a,b):
" Integral from ath data to bth data "
dim = b-a
if dim%2 != 0:
return "Error"
else:
x = data[a-1,0]
#xy = data[(b+a)/2-1,0]
y = data[b-1,0]
fx = data[a-1,0]
fy = data[b-1,0]
fxy = data[(b+a)/2-1,0]
return (y-x)/6*(fx+4*fxy+fy)```

### 复化求积公式/Composite Integration

#### 复化梯形公式

##### Composite Trapezium

```def com_trapezium(data):
dim = data.shape[0]
a = data[0,0]
b = data[dim-1,0]
h = (b-a)/(dim-1)
result = data[0,1]+data[dim-1,1]
for index in range(1,dim-1):
result = result + 2*data[index,1]
return result*h/2```

#### 复化Sinpson公式

##### Composite Sinpson

```def com_sinpson(data):
dim = data.shape[0]
a = data[0,0]
b = data[dim-1,0]
h = (b-a)/(dim-1)*2
result = data[0,1]+data[dim-1,1]
#times = (dim-1)/2
for index in range(1,dim-1):
if index%2 != 0:
result = result + 4*data[index,1]
else:
result = result + 2*data[index,1]
return result*h/6```

```data = np.array([[0,1],
[0.125,0.997398],
[0.25,0.989688],
[0.375,0.976727],
[0.5,0.958851],
[0.625,0.936156],
[0.75,0.908858],
[0.875,0.877193],
[1,0.841471]])
com_trapezium(data)
com_sinpson(data)```

``````0.94570081249999993
0.94609004166666677
``````

## 常微分方程数值解/Numerical Solution in Ordinary Differential Equation

### Euler法

#### Euler Method

Euler法使用折现逼近曲线，通过初值条件逐次迭代得出整条曲线

```import numpy as np
import matplotlib.pyplot as plt

""" Input:
Interval: [a,b]
Step: h
Differential: f(x,y)
Initial point: y(a) = y0
"""

def EulerMethod(a,b,h,f,y0):
if b>a:
step = int((b-a)/h+1)
x = np.linspace(a,b,step)
y = np.zeros_like(x)
y[0]=y0
for index in range(1,step):
y[index] = y[index-1] + h*f(x[index-1],y[index-1])
return y
else:
return "Error"```

### 改进的Euler法

#### Improved Euler Method

```def I_EulerMethod(a,b,h,f,y0):
if b>a:
step = int((b-a)/h+1)
x = np.linspace(a,b,step)
y = np.zeros_like(x)
ye = np.zeros_like(x)
y[0] = y0
ye[0] = y0
for index in range(1,step):
ye[index] = ye[index-1] + h*f(x[index-1],ye[index-1])
y[index] = y[index-1] + 0.5*h*(f(x[index-1],y[index-1])+f(x[index],ye[index]))
return y
else:
return "Error"```

```def Draw(a,b,h,f,y0,ft):
step = int((b-a)/h+1)
x = np.linspace(a,b,step)
fig, ax = plt.subplots(dpi=120)
ye = EulerMethod(a,b,h,f,y0)
yie = I_EulerMethod(a,b,h,f,y0)
yt = ft(x)
ax.plot(x,ye,marker='o',color='r',alpha=0.5,label='Euler Method')
ax.plot(x,yie,marker='x',color='b',alpha=0.5,label='Improved Euler Method')
ax.plot(x,yt,marker='v',color='g',alpha=0.5,label='Explicit Function')
ax.legend()
fig.show()
return 0```

```a = 0
b = 1
h = 0.1
f = lambda x,y: y-2*x/y
y0 = 1
ft = lambda x: np.sqrt(1+2*x)

EulerMethod(a,b,h,f,y0)
I_EulerMethod(a,b,h,f,y0)
Draw(a,b,h,f,y0,ft)```

``````array([ 1.        ,  1.1       ,  1.19181818,  1.27743783,  1.3582126 ,
1.43513292,  1.50896625,  1.58033824,  1.64978343,  1.71777935,
1.78477083])
array([ 1.        ,  1.09590909,  1.18438953,  1.26711005,  1.34524979,
1.41969469,  1.49114658,  1.56018901,  1.62733006,  1.69303203,
1.7577335 ])
``````
You can’t perform that action at this time.