求根的近似值
• 二分法，普通迭代，牛顿迭代三合一
"""
author:cxb
numberid:20210511085
time:20240319
"""
def bitsect(f, a, b, ep):
fa = f(a);
fb = f(b)
if fa * fb > 0:
index = 0;
return index;
k = 1
while abs(b - a) / 2 > ep:
x = (a + b) / 2;
fx = f(x)
if fx == 0:
index = 1;
return index, k, x
elif fx * fa < 0:
b = x;
fb = fx
else:
a = x;
fa = fx
k = k + 1
index = 1
return index, k, (a + b) / 2
def Newton(f, fdao, x, ep, it_max):
index = 0;
k = 1
while k < it_max:
x1 = x;
fv = f(x);
fdaov = fdao(x)
if abs(fdaov) < ep:
break
x = x1 - fv / fdaov
k = k + 1
if abs(x - x1) < ep:
index = 1
break
return index, k, x
def easyfind(f, x, ep):
flag = 0;
k = 1;
last = x
#改动
x = (3*x + 1) ** (1 / 3)
#改动
while abs(x - last) > ep:
last = x;
#改动
x = (3*x + 1) ** (1 / 3)
#改动
k = k + 1
if (abs(x - last) <= ep):
flag = 1
return flag, k, x
#改动 把自己的函数写进去
f = lambda x: x ** 3 +0*x ** 2-3*x - 1
fdao = lambda x: 3 * x ** 2 + 0* x - 3
#改动 f不改，前两个数是区间范围，第三个数是二分次数，最后是精度
jieguo_b = bitsect(f, 1, 2, 2,0.001)
if jieguo_b[0] == 1:
print('successed to binary search in {} times,the x* is {}'.format(jieguo_b[1], jieguo_b[2]))
else:
print('failed to binary search')
#改动 x0 精度 次数
jieguo_N = Newton(f, fdao, 2,1, 3)
if jieguo_N[0] == 1:
print('successed to Newton-I in {} times,the x* is {}'.format(jieguo_N[1], jieguo_N[2]))
else:
print('failed to Newton-I')
#改动 1应该是题目会给的初值x0，中间是次数，题目后面会给精度（运行不了就降低精度）
easyjieguo = easyfind(f, 2,3, 0.01)
if easyjieguo[0] == 1:
print('successed to eayfind in {} times,the x* is {}'.format(easyjieguo[1], easyjieguo[2]))
else:
print('failed to easyfind search')
高斯消元
"""
Created on Tue Mar 26 15:12:22 2024
@author: cxb
"""
import numpy as np
exp = 1e-7
def Gauss(A, b):
n = len(b)
index = 1
x = np.zeros(n)
for k in range(n):
# '选主元'
a_max = 0
for i in range(k, n):
if abs(A[i][k]) > a_max:
a_max = abs(A[i][k])
r = i
if a_max < exp:
index = 0
return
# '交换两行'
if r > k:
for j in range(k, n):
z = A[k][j]
A[k][j] = A[r][j]
A[r][j] = z
z = b[k]
b[k] = b[r]
b[r] = z
# 消元计算
for i in range(k+1, n):
m = A[i][k]/A[k][k]
for j in range(k+1, n):
A[i][j] -= m*A[k][j]
b[i] -= m*b[k]
# 回代过程
if abs(A[n-1][n-1]) < exp:
index = 0
return
for k in range(n-1, -1, -1):
for j in range(k+1, n):
b[k] -= A[k][j]*x[j]
x[k] = b[k]/A[k][k]
return index, x
#改动
A = [[2, 1, 0], [1, 3, 1], [0, 1, 2]]
b = [1, -1, 1]
print(“Gauss:”,end=‘’)
Print(Gauss(A,b,1e-7))
雅可比迭代
"""
Created on Tue Apr 2 15:13:50 2024
@Aurhor :Xiangbo Cheng
"""
import numpy as np
def Jacobi(A, b, x0, it_max, ep):
D = np.diag(np.diag(A));
U = -np.triu(A, 1);
L = -np.tril(A, -1)
B = np.dot(np.linalg.inv(D), (L + U))
f = np.dot(np.linalg.inv(D), b)
x = np.dot(B, x0) + f
k = 1;
index = 1
while np.linalg.norm(x - x0) >= ep:
x0 = x
x = np.dot(B, x0) + f
k = k + 1
if k > it_max:
index = 0;
break
return k, index, x
#改动 it_max 是次数，一次只能手算，两次改为1，三次改为2
A = [[2, 1, 0], [1, 3, 1], [0, 1, 2]]
b = [1, -1, 1]
x0 = [0, 0, 0]
print(“Jacobi:”,end=‘’)
print(Jacobi(A, b, x0, 2, 1e-6))
实验数据
A = [[15, 3, 5, 3], [2, 15, 6, 5], [3, 6, 18, 4], [2, 5, 3, 13]]
b = [12, 8, 20, 5]
x0 = [0, 0, 0, 0]
print(Jacobi(A, b, x0, 100, 1e-6))
高斯-赛德尔迭代
"""
Created on Tue Apr 2 15:13:50 2024
@Aurhor :Xiangbo Cheng
"""
import numpy as np
def Gauss_Seidel(A, b, x0, it_max, ep):
D = np.diag(np.diag(A));
U = -np.triu(A, 1);
L = -np.tril(A, -1)
B = np.dot(np.linalg.inv(D - L), (U))
f = np.dot(np.linalg.inv(D - L), b)
x = np.dot(B, x0) + f
k = 1;
index = 1
while np.linalg.norm(x - x0) >= ep:
x0 = x
x = np.dot(B, x0) + f
k = k + 1
if k > it_max:
index = 0;
break
return k, index, x
#改动
A = [[2, 1, 0], [1, 3, 1], [0, 1, 2]]
b = [1, -1, 1]
x0 = [0, 0, 0]
print(“Gauss_Seidel:”,end=‘’)
print(Gauss_Seidel(A, b, x0, 2, 0.001))
实验数据
A = [[15, 3, 5, 3], [2, 15, 6, 5], [3, 6, 18, 4], [2, 5, 3, 13]]
b = [12, 8, 20, 5]
x0 = [0, 0, 0, 0]
print(Gauss_Seidel(A, b, x0, 100, 1e-6))
拉格朗日插值法
# -*- coding: utf-8 -*-
"""
Created on Sun Apr 7 15:06:43 2024
@author: JiangHao
"""
from math import sqrt
import pylab as pl
def LagLangRi(x, y, xi):
n = len(x)
p = [1] * n
yi = 0
for k in range(n):
t = [1] * n
for j in range(n):
if j != k:
if abs(x[k] - x[j]) < 0.0000001:
return 0
t[j] = (xi - x[j]) / (x[k] - x[j])
p[k] = p[k] * t[j]
yi = yi + y[k] * p[k]
return yi
'''
x=[1,4,9,16]
y=[1,2,3,4]
xxx=[];yyy=[];yyy1=[]
for i in range(10,160):
xxx.append(i/10)
yyy.append(LagLangRi(x,y,i/10))
yyy1.append(sqrt(i/10))
pl.plot(xxx,yyy,label="LagLangRi")
pl.plot(xxx,yyy1,label="real")
pl.legend()
'''
x = [0] * 6
y = [0] * 6
for i in range(6):
x[i] = -5 + i * 2
y[i] = 1 / (1 + x[i] ** 2)
x10 = [0] * 11
y10 = [0] * 11
for i in range(11):
x10[i] = -5 + i * 1
y10[i] = 1 / (1 + x10[i] ** 2)
xxx = [0] * 101
yyy = [0] * 101
yyy1 = [0] * 101
yyy2 = [0] * 101
for i in range(101):
xxx[i] = -5 + 0.1 * i
yyy[i] = LagLangRi(x, y, xxx[i])
yyy1[i] = LagLangRi(x10, y10, xxx[i])
yyy2[i] = 1 / (1 + xxx[i] ** 2)
pl.plot(xxx, yyy, label="LagLangRi-5")
pl.plot(xxx, yyy1, label="LagLangRi-10")
pl.plot(xxx, yyy2, label="real")
pl.legend()
分段线性插值法
# -*- coding: utf-8 -*-
"""
Created on Tue Apr 9 15:13:50 2024
@author: cxb
"""
import pylab as pl
from math import sqrt
def Line_int(x, y, xi):
n = len(x)
for k in range(n - 1):
if (abs(x[k] - x[k + 1])) < 1e-7:
return 0
if (x[k] <= xi <= x[k + 1]):
y1 = (xi - x[k + 1]) / (x[k] - x[k + 1]) * y[k] + (xi - x[k]) / (x[k + 1] - x[k]) * y[k + 1]
return y1
x = [1, 4, 9, 16]
y = [1, 2, 3, 4]
xxx = []
yyy1 = []
yyy2 = []
for i in range(10, 160):
xxx.append(i / 10)
yyy1.append(sqrt(i / 10))
yyy2.append(Line_int(x, y, i / 10))
pl.plot(xxx, yyy1, label="y=sqrt(x)", linestyle='--')
pl.plot(xxx, yyy2, label="line-int")
pl.legend()
两种插值法的比对
# -*- coding: utf-8 -*-
"""
Created on Sun Apr 7 15:06:43 2024
@author: Administrator
"""
from math import sqrt
import pylab as pl
def Line_int(x, y, xi):
n = len(x)
for k in range(n - 1):
if (abs(x[k] - x[k + 1])) < 1e-7:
return 0
if (x[k] <= xi <= x[k + 1]):
y1 = (xi - x[k + 1]) / (x[k] - x[k + 1]) * y[k] + (xi - x[k]) / (x[k + 1] - x[k]) * y[k + 1]
return y1
def LagLangRi(x, y, xi):
n = len(x)
p = [1] * n
yi = 0
for k in range(n):
t = [1] * n
for j in range(n):
if j != k:
if abs(x[k] - x[j]) < 1e-7:
return 0
t[j] = (xi - x[j]) / (x[k] - x[j])
p[k] = p[k] * t[j]
yi = yi + y[k] * p[k]
return yi
'''
x=[1,4,9,16]
y=[1,2,3,4]
xxx=[];yyy=[];yyy1=[]
for i in range(10,160):
xxx.append(i/10)
yyy.append(LagLangRi(x,y,i/10))
yyy1.append(sqrt(i/10))
pl.plot(xxx,yyy,label="LagLangRi")
pl.plot(xxx,yyy1,label="real")
pl.legend()
'''
x = [0] * 6
y = [0] * 6
for i in range(6):
x[i] = -5 + i * 2
y[i] = 1 / (1 + x[i] ** 2)
x10 = [0] * 11
y10 = [0] * 11
for i in range(11):
x10[i] = -5 + i * 1
y10[i] = 1 / (1 + x10[i] ** 2)
xxx = [0] * 101
yyy = [0] * 101
yyy1 = [0] * 101
yyy2 = [0] * 101
yyy3 = [0] * 101
for i in range(101):
xxx[i] = -5 + 0.1 * i
yyy[i] = LagLangRi(x, y, xxx[i])
yyy1[i] = LagLangRi(x10, y10, xxx[i])
yyy2[i] = 1 / (1 + xxx[i] ** 2)
yyy3[i] = Line_int(x10, y10, xxx[i])
pl.plot(xxx, yyy, label="LagLangRi-5")
pl.plot(xxx, yyy1, label="LagLangRi-10")
pl.plot(xxx, yyy2, label="real")
pl.plot(xxx, yyy3, label="Line-int")
pl.legend()
三次样条插值法
# -*- coding: utf-8 -*-
import pylab as pl
import numpy as np
def cubic_spline2(x, y, ydot, xi):
n = len(x);
mm = len(xi);
yi = np.zeros(mm)
eps = 1e-6;
if ydot == []:
ydot = [(y[1] - y[0]) / (x[1] - x[0]), (y[n - 1] - y[n - 2]) / (x[n - 1] - x[n - 2])];
h = np.zeros(n);
lamda = np.ones(n);
mu = np.ones(n - 2);
m = np.zeros(n - 2);
d = np.zeros(n - 2);
for k in range(1, n):
h[k] = x[k] - x[k - 1];
if abs(h[k]) < eps:
return 0;
for k in range(1, n - 1):
lamda[k - 1] = h[k + 1] / (h[k] + h[k + 1]);
mu[k - 1] = 1 - lamda[k - 1];
d[k - 1] = 3 * (mu[k - 1] * (y[k + 1] - y[k]) / h[k + 1] + lamda[k - 1] * (y[k] - y[k - 1]) / h[k]);
d[0] = d[0] - lamda[0] * ydot[0];
d[n - 3] = d[n - 3] - mu[n - 3] * ydot[1];
A = np.diag(2 * np.ones(n - 2))
for i in range(0, n - 3):
A[i, i + 1] = mu[i];
A[i + 1, i] = lamda[i + 1];
m = np.dot(np.linalg.inv(A), d)
m = np.insert(m, 0, ydot[0])
m = np.insert(m, n - 1, ydot[1])
for i in range(0, mm):
for k in range(1, n):
if x[k - 1] <= xi[i] <= x[k]:
yi[i] = y[k - 1] / h[k] ** 3 * (xi[i] - x[k]) ** 2 * (h[k] + 2 * (xi[i] - x[k - 1])) + y[k] / h[
k] ** 3 * (xi[i] - x[k - 1]) ** 2 * (h[k] + 2 * (x[k] - xi[i])) \
+ m[k - 1] / h[k] ** 2 * (xi[i] - x[k - 1]) * (xi[i] - x[k]) ** 2 + m[k] / h[k] ** 2 * (
xi[i] - x[k]) * (xi[i] - x[k - 1]) ** 2
break
return yi
x10 = [0] * 11;
y10 = [0] * 11
for i in range(11):
x10[i] = -5 + i * 1
y10[i] = 1 / (1 + x10[i] ** 2)
xxx = [0] * 101;
yyy1 = [0] * 101;
yyy2 = [0] * 101
for i in range(101):
xxx[i] = -5 + 0.1 * i
yyy1[i] = 1 / (1 + xxx[i] ** 2)
yyy2 = cubic_spline2(x10, y10, [], xxx)
pl.plot(xxx, yyy2, label="spline")
# pl.plot(xxx,yyy1,label="real")
pl.legend()
最小二乘
# -*- coding: utf-8 -*-
"""
Created on Thu Apr 11 15:11:33 2024
@author: Administrator
"""
import numpy as np
import pylab as pl
def phi_k(x, k):
y = x ** k
return y
def squar_least(x, y, w, n):
G = np.zeros((n + 1, n + 1))
d = np.zeros((n + 1))
for i in range(n + 1):
for j in range(n + 1):
G[i, j] = sum((w * phi_k(x, i)) * phi_k(x, j))
for i in range(n + 1):
d[i] = sum((w * phi_k(x, i)) * y)
a = np.dot(np.linalg.inv(G), d)
return a
x = np.array([0, 1, 4.04, 9, 16.1, 25, 36, 49.2, 64.2])
y = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8])
w = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1])
a = squar_least(x, y, w, 4)
yuce_y = a[0] * 1 + a[1] * 0.6 + a[2] * 0.6 * 0.6 + a[3] * 0.6 ** 3 + a[4] * 0.6 ** 4
print(yuce_y)
xxx = []
yyy = []
for i in range(0, 120):
x1 = i / 100
xxx.append(x1)
yyy.append(a[0] * 1 + a[1] * x1 + a[2] * x1 ** 2 + a[3] * x1 ** 3 + a[4] * x1 ** 4)
pl.scatter(x, y)
pl.plot(xxx, yyy)
复合梯形求积公式
# -*- coding: utf-8 -*-
"""
Created on Tue Apr 16 15:20:16 2024
@author: ChengXiangbo
"""
import numpy as np
import math
def T_quad(x, y):
n = len(x)
m = len(y)
if n != m:
return 0
h = (x[n - 1] - x[0]) / (n - 1)
sum_T = 0
for i in range(n - 1):
sum_T = sum_T + h / 2 * (y[i] + y[i + 1])
return sum_T
x = np.arange(0, 1 + 1 / 8, 0.125)
y = []
for i in x:
if i == 0:
y.append(1)
else:
y.append(math.sin(i) / i)
quadT = T_quad(x, y)
print(quadT)
高斯复合求积公式
# -*- coding: utf-8 -*-
"""
Created on Tue Apr 16 15:20:16 2024
@author: Administrator
"""
from math import exp
from math import sqrt
from math import log
f = lambda x: 1 / x
def G_quad(f, a, b, N):
h = (b - 1) / N;
t = [-sqrt(3 / 5), 0, sqrt(3 / 5)]
A = [5 / 9, 8 / 9, 5 / 9]
x_fromt = [0, 0, 0]
Sum_G = 0
for k in range(N):
for i in range(3):
x_fromt[i] = h / 2 * t[i] + a + (k + 1 / 2) * h
Sum_G = Sum_G + h / 2 * (A[0] * f(x_fromt[0]) + A[1] * f(x_fromt[1]) + A[2] * f(x_fromt[2]))
return Sum_G
quadG = G_quad(f, 1, 5, 15)
print(quadG)
print(log(5) - quadG)
辛普森复合求积公式
# -*- coding: utf-8 -*-
"""
Created on Tue Apr 16 15:20:16 2024
@author: Xiangbo Cheng
"""
import numpy as np
import math
def fun(x):
if x == 0:
return 1
return math.sin(x) / x
def xps(a, b, n):
h = (b - a) / n
x = a
s = fun(x) - fun(b)
for i in range(1, n + 1):
x = x + h / 2
s = s + 4 * fun(x)
x = x + h / 2
s = s + 2 * fun(x)
result = (h / 6) * s
return result
a = 0
b = 1
p = 4
print(xps(a, b, p))
常微分方程的数值解法
# -*- coding: utf-8 -*-
"""
Created on Tue Apr 23 15:03:04 2024
@author: cxb
"""
from math import *
import numpy as np
import pylab as pl
f = lambda x, y: y - 2 * x / y
def Euler(f, x0, y0, h, N):
x = np.zeros(N + 1)
y = np.zeros(N + 1)
x[0] = x0
y[0] = y0
for i in range(N):
x[i + 1] = x[i] + h
y[i + 1] = y[i] + h * f(x[i], y[i])
return x, y
def Rounge_kutta4(f, x0, y0, h, N):
x = np.zeros(N + 1)
y = np.zeros(N + 1)
x[0] = x0
y[0] = y0
for i in range(N):
x[i + 1] = x[i] + h
k1 = f(x[i], y[i])
k2 = f(x[i] + h / 2, y[i] + h / 2 * k1)
k3 = f(x[i] + h / 2, y[i] + h / 2 * k2)
k4 = f(x[i] + h, y[i] + h * k3)
y[i + 1] = y[i] + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4)
return x, y
def Gai_Euler(f, x0, y0, h, N):
x = np.zeros(N + 1)
y = np.zeros(N + 1)
x[0] = x0
y[0] = y0
for i in range(N):
x[i + 1] = x[i] + h
y[i + 1] = y[i] + h * f(x[i], y[i])
y[i + 1] = y[i] + h * f(x[i + 1], y[i + 1])
return x, y
[xx_r, yy_r] = Euler(f, 0, 1, 0.1, 10)
[xx_Gai, yy_Gai] = Gai_Euler(f, 0, 1, 0.1, 10)
[xx_rk, yy_rk] = Rounge_kutta4(f, 0, 1, 0.1, 10)
xxx = np.arange(0, 1 + 0.01, 0.01)
yyy = []
for i in xxx:
yyy.append(sqrt(1 + 2 * i))
pl.plot(xx_r, yy_r, label="Euler", color="red")
pl.plot(xx_Gai, yy_Gai, label="Gai_Euler", color='blue')
pl.plot(xx_rk, yy_rk, label="Rounge_kutta4", color="green")
pl.plot(xxx, yyy, linestyle=':', label="y=sqrt(2x+1)", color='yellow')
pl.legend()
1. 用这9个点作8次多项式插值
# -*- coding: utf-8 -*-
"""
Created on Sun Apr 7 15:06:43 2024
@author: JiangHao
"""
from math import sqrt
import numpy as np
import pylab as pl
def LagLangRi(x, y, xi):
n = len(x)
p = [1] * n
yi = 0
for k in range(n):
t = [1] * n
for j in range(n):
if j != k:
if abs(x[k] - x[j]) < 0.0000001:
return 0
t[j] = (xi - x[j]) / (x[k] - x[j])
p[k] = p[k] * t[j]
yi = yi + y[k] * p[k]
return yi
x = np.array([0, 1, 4.04, 9, 16.1, 25, 36, 49.2, 64.2])
y = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8])
xxx=[];yyy=[];yyy1=[]
for i in range(0,640):
xxx.append(i/10)
yyy.append(LagLangRi(x,y,i/10))
yyy1.append(sqrt(i/10))
pl.plot(xxx,yyy,label="LagLangRi")
pl.plot(xxx,yyy1,label="real")
pl.legend()
2.进行分段线性插值
# -*- coding: utf-8 -*-
"""
Created on Tue Apr 9 15:13:50 2024
@author: cxb
"""
from math import sqrt
import numpy as np
import pylab as pl
def Line_int(x, y, xi):
n = len(x)
for k in range(n - 1):
if (abs(x[k] - x[k + 1])) < 1e-7:
return 0
if (x[k] <= xi <= x[k + 1]):
y1 = (xi - x[k + 1]) / (x[k] - x[k + 1]) * y[k] + (xi - x[k]) / (x[k + 1] - x[k]) * y[k + 1]
return y1
x = np.array([0, 1, 4.04, 9, 16.1, 25, 36, 49.2, 64.2])
y = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8])
xxx = []
yyy1 = []
yyy2 = []
for i in range(0, 642):
xxx.append(i / 10)
yyy1.append(sqrt(i / 10))
yyy2.append(Line_int(x, y, i / 10))
pl.plot(xxx, yyy1, label="y=sqrt(x)", linestyle='--')
pl.plot(xxx, yyy2, label="line-int")
pl.legend()
3. 进行三次样条插值；
# -*- coding: utf-8 -*-
from math import sqrt
import numpy as np
import pylab as pl
def cubic_spline2(x, y, ydot, xi):
n = len(x);
mm = len(xi);
yi = np.zeros(mm)
eps = 1e-6;
if ydot == []:
ydot = [(y[1] - y[0]) / (x[1] - x[0]), (y[n - 1] - y[n - 2]) / (x[n - 1] - x[n - 2])];
h = np.zeros(n);
lamda = np.ones(n);
mu = np.ones(n - 2);
m = np.zeros(n - 2);
d = np.zeros(n - 2);
for k in range(1, n):
h[k] = x[k] - x[k - 1];
if abs(h[k]) < eps:
return 0;
for k in range(1, n - 1):
lamda[k - 1] = h[k + 1] / (h[k] + h[k + 1]);
mu[k - 1] = 1 - lamda[k - 1];
d[k - 1] = 3 * (mu[k - 1] * (y[k + 1] - y[k]) / h[k + 1] + lamda[k - 1] * (y[k] - y[k - 1]) / h[k]);
d[0] = d[0] - lamda[0] * ydot[0];
d[n - 3] = d[n - 3] - mu[n - 3] * ydot[1];
A = np.diag(2 * np.ones(n - 2))
for i in range(0, n - 3):
A[i, i + 1] = mu[i];
A[i + 1, i] = lamda[i + 1];
m = np.dot(np.linalg.inv(A), d)
m = np.insert(m, 0, ydot[0])
m = np.insert(m, n - 1, ydot[1])
for i in range(0, mm):
for k in range(1, n):
if x[k - 1] <= xi[i] <= x[k]:
yi[i] = y[k - 1] / h[k] ** 3 * (xi[i] - x[k]) ** 2 * (h[k] + 2 * (xi[i] - x[k - 1])) + y[k] / h[
k] ** 3 * (xi[i] - x[k - 1]) ** 2 * (h[k] + 2 * (x[k] - xi[i])) \
+ m[k - 1] / h[k] ** 2 * (xi[i] - x[k - 1]) * (xi[i] - x[k]) ** 2 + m[k] / h[k] ** 2 * (
xi[i] - x[k]) * (xi[i] - x[k - 1]) ** 2
break
return yi
x = np.array([0, 1, 4.04, 9, 16.1, 25, 36, 49.2, 64.2])
y = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8])
xxx=[];yyy=[];yyy1=[]
for i in range(0,640):
xxx.append(i/10)
yyy.append(cubic_spline2(x,y,[],xxx))
yyy1.append(sqrt(i/10))
pl.plot(xxx, yyy1, label="spline")
pl.plot(xxx,yyy1,label="real")
pl.legend()
4.用最小二乘拟合方法进行二次多项式拟合
# -*- coding: utf-8 -*-
"""
Created on Thu Apr 11 15:11:33 2024
@author: Administrator
"""
import numpy as np
import pylab as pl
def phi_k(x, k):
y = x ** k
return y
def squar_least(x, y, w, n):
G = np.zeros((n + 1, n + 1))
d = np.zeros((n + 1))
for i in range(n + 1):
for j in range(n + 1):
G[i, j] = sum((w * phi_k(x, i)) * phi_k(x, j))
for i in range(n + 1):
d[i] = sum((w * phi_k(x, i)) * y)
a = np.dot(np.linalg.inv(G), d)
return a
x = np.array([0, 1, 4.04, 9, 16.1, 25, 36, 49.2, 64.2])
y = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8])
w = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1])
a = squar_least(x, y, w, 4)
yuce_y = a[0] * 1 + a[1] * 0.6 + a[2] * 0.6 * 0.6 + a[3] * 0.6 ** 3 + a[4] * 0.6 ** 4
print(yuce_y)
xxx = []
yyy = []
for i in range(0, 6420):
x1 = i / 100
xxx.append(x1)
yyy.append(a[0] * 1 + a[1] * x1 + a[2] * x1 ** 2 + a[3] * x1 ** 3 + a[4] * x1 ** 4)
pl.scatter(x, y)
pl.plot(xxx, yyy)
实验四
1. 复合梯形
# -*- coding: utf-8 -*-
"""
Created on Tue Apr 16 15:20:16 2024
@author: ChengXiangbo
"""
import numpy as np
import math
def T_quad(x, y):
n = len(x)
m = len(y)
if n != m:
return 0
h = (x[n - 1] - x[0]) / (n - 1)
sum_T = 0
for i in range(n - 1):
sum_T = sum_T + h / 2 * (y[i] + y[i + 1])
return sum_T
x = np.arange(1, 9+0.1, 0.1)
y = []
for i in x:
y.append(1/i)
quadT = T_quad(x, y)
print(quadT)
2. 辛普森
# -*- coding: utf-8 -*-
"""
Created on Tue Apr 16 15:20:16 2024
@author: Xiangbo Cheng
"""
import numpy as np
import math
def fun(x):
if x == 0:
return 1
return 1 / x
def xps(a, b, n):
h = (b - a) / n
x = a
s = fun(x) - fun(b)
for i in range(1, n + 1):
x = x + h / 2
s = s + 4 * fun(x)
x = x + h / 2
s = s + 2 * fun(x)
result = (h / 6) * s
return result
a = 1
b = 9
p = 4
print(xps(a, b, p))
4. 高斯复合
# -*- coding: utf-8 -*-
"""
Created on Tue Apr 16 15:20:16 2024
@author: Administrator
"""
from math import exp
from math import sqrt
from math import log
f = lambda x: 1 / x
def G_quad(f, a, b, N):
h = (b - 1) / N;
t = [-sqrt(3 / 5), 0, sqrt(3 / 5)]
A = [5 / 9, 8 / 9, 5 / 9]
x_fromt = [0, 0, 0]
Sum_G = 0
for k in range(N):
for i in range(3):
x_fromt[i] = h / 2 * t[i] + a + (k + 1 / 2) * h
Sum_G = Sum_G + h / 2 * (A[0] * f(x_fromt[0]) + A[1] * f(x_fromt[1]) + A[2] * f(x_fromt[2]))
return Sum_G
quadG = G_quad(f, 1, 9, 20)
print(quadG)
print(log(9) - quadG)
实验五
# -*- coding: utf-8 -*-
"""
Created on Tue Apr 23 15:03:04 2024
@author: Administrator
"""
from math import *
import numpy as np
import pylab as pl
f = lambda x, y: -50*y+50*x*x+2*x
def Euler(f, x0, y0, h, N):
x = np.zeros(N + 1)
y = np.zeros(N + 1)
x[0] = x0
y[0] = y0
for i in range(N):
x[i + 1] = x[i] + h
y[i + 1] = y[i] + h * f(x[i], y[i])
return x, y
def Rounge_kutta4(f, x0, y0, h, N):
x = np.zeros(N + 1)
y = np.zeros(N + 1)
x[0] = x0
y[0] = y0
for i in range(N):
x[i + 1] = x[i] + h
k1 = f(x[i], y[i])
k2 = f(x[i] + h / 2, y[i] + h / 2 * k1)
k3 = f(x[i] + h / 2, y[i] + h / 2 * k2)
k4 = f(x[i] + h, y[i] + h * k3)
y[i + 1] = y[i] + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4)
return x, y
def Gai_Euler(f, x0, y0, h, N):
x = np.zeros(N + 1)
y = np.zeros(N + 1)
x[0] = x0
y[0] = y0
for i in range(N):
x[i + 1] = x[i] + h
y[i + 1] = y[i] + h * f(x[i], y[i])
y[i + 1] = y[i] + h * f(x[i + 1], y[i + 1])
return x, y
[xx_r, yy_r] = Euler(f, 0, 1/2, 0.01, 100)
[xx_Gai, yy_Gai] = Gai_Euler(f, 0, 1/2, 0.01, 100)
[xx_rk, yy_rk] = Rounge_kutta4(f, 0, 1/2, 0.01, 100)
xxx = np.arange(0, 1 + 0.01, 0.01)
yyy = []
for i in xxx:
yyy.append(e**(-50*i)/3+i**2)
pl.plot(xx_r, yy_r, label="Euler", color="red")
pl.plot(xx_Gai, yy_Gai, label="Gai_Euler", color='blue')
pl.plot(xx_rk, yy_rk, label="Rounge_kutta4", color="green")
pl.plot(xxx, yyy, linestyle=':', label="y(x)=e-50x/3+x2", color='yellow')
pl.legend()
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
