#二分法 a和b是区间，ep是要求的精度
def bitsect(f, a, b, it_max,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
if k > it_max:
index = 0;
break
index = 1
return index, k, (a + b) / 2
# def easyfind(f, x, ep):
# flag = 0;
# k = 1;
# last = x
# x = (x + 1) ** (1 / 3)
# while abs(x - last) > ep:
# last = x;
# x = (x + 1) ** (1 / 3)
# k = k + 1
# if (abs(x - last) <= ep):
# flag = 1
# return flag, k, x
# #简单迭代法 x是题目给的初值,ep是精度
def easyfind(f, x, it_max,ep):
flag = 0;
k = 1;
last = x
# 改动 x是转换的f(x),例如f=x**3-x-1,x=(x + 1) ** (1 / 3)
x = (3*x + 1) ** (1 / 3)
# 改动
while abs(x - last) > ep:
last = x;
# 改动 x是转换的f(x),例如f=x**3-x-1,x=(x + 1) ** (1 / 3)
x = (3*x + 1) ** (1 / 3)
# 改动
k = k + 1
if k > it_max:
break
if (abs(x - last) <= ep):
flag = 1
return flag, k, x
#牛顿迭代法 x是取的初值可以0.5、1、1.5，ep是精度要求，it max是最多迭代多少次一般取100
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
# 改动
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]-1, jieguo_b[2]))
else:
print('failed to binary search')
# 改动 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')
# 改动 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. 二分法：
# ○ 优点：二分法是一种简单而稳定的迭代求根方法。它适用于连续、单调函数的求根问题，并且对于函数值异号的根有较好的收敛性。二分法的收敛速度相对较慢，但每次迭代可以保证区间缩小一半，因此可以保证最终达到指定的精度要求。
# ○ 缺点：二分法的收敛速度比较慢，尤其是对于复杂非线性函数，可能需要较多的迭代次数才能达到较高的精度要求。此外，二分法只能应用于单变量函数的求根问题，对于多变量函数需要进行相应的扩展。
# 2. 普通迭代法（也称为不动点迭代法）：
# ○ 优点：普通迭代法是一种简单而直观的迭代求根方法。它通过将方程转化为不动点形式，通过不断迭代逼近不动点从而求解根。普通迭代法通常适用于函数局部收敛的情况，并且可以通过简单的迭代公式进行实现。
# ○ 缺点：普通迭代法的收敛性和收敛速度高度依赖于迭代函数的选择。对于某些函数，可能出现发散的情况，或者收敛速度较慢。此外，普通迭代法的收敛性还与初始迭代点的选择有关。
# 3. 牛顿迭代法：
# ○ 优点：牛顿迭代法是一种快速收敛的迭代求根方法，尤其适用于函数局部收敛的情况。它通过利用函数的一阶导数信息，通过迭代逼近函数的根。牛顿迭代法的收敛速度通常是二阶的，因此它可以迅速逼近根，并且对于复杂非线性函数也有较好的效果。
# ○ 缺点：牛顿迭代法的收敛性和收敛速度高度依赖于初始迭代点的选择。在某些情况下，初始点选择不当可能导致迭代发散或陷入局部极值。此外，牛顿迭代法要求函数的一阶导数存在且不为零，因此不适用于某些函数或者在不好求导的情况下。
# 综上所述，二分法适用于连续函数中根的求解，稳定但收敛速度较慢；普通迭代法简单直观，但收敛性和收敛速度依赖于迭代函数的选择；牛顿迭代法收敛速度快，适用于函数局部收敛的情况，但初始点选择和导数的存在性对其影响较大。在实际应用中，根据函数特点选择适当的求根方法是很重要的。
下是二///
import numpy as np
#高斯列主元消去法 A，b是矩阵，exp是精度
def Gauss(A, b,exp):
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
# 改动 1e-7就是精度，一般不用改，或者看题目要求
A = [[2,1,0], [1,3,1], [0,1,2]]
b = [1,-1,1]
print("Gauss:",end='')
print(Gauss(A, b,1e-7))
# 改动
#雅可比迭代法A，b是矩阵，xo是题目给的初值，it max是最多迭代数，一般就是100，ep是精度例如0.0001
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...题目没要求就100 1e-6就是精度，一般不用改，或者看题目要求
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，b是矩阵，xo是题目给的初值，it max是最多迭代数，一般就是100，如需改动和上面一样 ep是精度例如0.00001
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))
# 改动
#高斯消元法、雅可比迭代法和高斯赛德尔迭代法是用来求解线性方程组的数值方法之一。它们都有各自的优点和限制，适用于不同的问题和情景。下面是这三种方法的简要介绍和优劣比较：
# 1. 高斯消元法：
# ○ 优点：高斯消元法是一种直接解法，能够精确求解线性方程组。它通过变换矩阵来消去未知元的系数，将方程组转化为上三角矩阵或对角矩阵，最后通过回代得到解。适用于小规模、稀疏的线性方程组求解。
# ○ 缺点：高斯消元法的计算复杂度较高，时间复杂度约为O(n^3)，对于大规模的线性方程组求解效率较低。当矩阵存在主元为零或接近零的情况时，可能出现除以零或数值精度损失等问题。
# 2. 雅可比迭代法：
# ○ 优点：雅可比迭代法是一种迭代解法，适用于大规模的稀疏线性方程组求解。它将线性方程组转化为迭代形式，通过不断更新未知元的估计值来逼近真实解。每个迭代步骤都可以在并行计算中进行，适用于分布式环境。
# ○ 缺点：雅可比迭代法通常需要进行大量的迭代才能达到较高的精度，收敛速度较慢。对于具有严格对角占优的线性方程组，收敛速度较快，但对于其他类型的方程组，可能需要更多的迭代次数。
# 3. 高斯赛德尔迭代法：
# ○ 优点：高斯赛德尔迭代法是雅可比迭代法的改进版本，它在每次迭代中利用更新后的未知元来提高收敛速度。相比于雅可比迭代法，高斯赛德尔迭代法通常能够更快地收敛到精确解。
# ○ 缺点：高斯赛德尔迭代法在进行更新时，需要使用已更新的未知元，因此无法进行并行计算。此外，对于某些线性方程组，高斯赛德尔迭代法可能发散或不收敛。
# 综上所述，高斯消元法适用于小规模、稀疏的线性方程组求解；雅可比迭代法适用于大规模的稀疏线性方程组求解，特别是具有严格对角占优性质的方程组；而高斯赛德尔迭代法在迭代收敛速度方面优于雅可比迭代法，但不适用于并行计算。选择适当的方法应根据具体问题的规模、稀疏性和收敛要求来决定。
//欧拉//
import numpy as np
import matplotlib.pyplot as plt
# 定义微分方程 dy/dt = -y
def dydt(y):
return -y
# 欧拉法
def euler(dydt, y0, t0, t_end, h):
t = np.arange(t0, t_end, h)
y = np.zeros_like(t)
y[0] = y0
for i in range(1, len(t)):
y[i] = y[i-1] + h * dydt(y[i-1])
return t, y
# 改进的欧拉法
def improved_euler(dydt, y0, t0, t_end, h):
t = np.arange(t0, t_end, h)
y = np.zeros_like(t)
y[0] = y0
for i in range(1, len(t)):
y_next = y[i-1] + h * dydt(y[i-1])
y[i] = y[i-1] + h * (dydt(y[i-1]) + dydt(y_next)) / 2
return t, y
# 四阶Runge-Kutta法
def runge_kutta(dydt, y0, t0, t_end, h):
t = np.arange(t0, t_end, h)
y = np.zeros_like(t)
y[0] = y0
for i in range(1, len(t)):
k1 = dydt(y[i-1])
k2 = dydt(y[i-1] + h/2 * k1)
k3 = dydt(y[i-1] + h/2 * k2)
k4 = dydt(y[i-1] + h * k3)
y[i] = y[i-1] + h/6 * (k1 + 2*k2 + 2*k3 + k4)
return t, y
# 初始条件和参数
y0 = 1
t0 = 0
t_end = 1
# 步长
hs = [0.1, 0.025, 0.01]
# 对每种方法和每个步长进行计算
for h in hs:
# 欧拉法
t_euler, y_euler = euler(dydt, y0, t0, t_end, h)
# 改进的欧拉法
t_improved, y_improved = improved_euler(dydt, y0, t0, t_end, h)
# Runge-Kutta法
t_rk, y_rk = runge_kutta(dydt, y0, t0, t_end, h)
# 绘制结果
plt.figure(figsize=(10, 6))
plt.plot(t_euler, np.exp(-t_euler), label='Exact')
plt.plot(t_euler, y_euler, 'ro-', label='Euler, h=' + str(h))
plt.plot(t_improved, y_improved, 'gx-', label='Improved Euler, h=' + str(h))
plt.plot(t_rk, y_rk, 'bv-', label='RK4, h=' + str(h))
plt.legend()
plt.title('Comparison of Numerical Methods for ODEs')
plt.xlabel('Time')
plt.ylabel('Solution')
plt.grid(True)
plt.show()
。插值///
def f(x):
# 定义被积函数，例如 f(x) = x^2
return x**2
def composite_trapezoidal_rule(a, b, n):
# a 和 b 是积分区间的端点，n 是小区间的数量
h = (b - a) / n # 计算步长
sum = (f(a) + f(b)) / 2 # 计算端点的函数值的平均
for i in range(1, n):
sum += f(a + i * h) # 累加每个小区间的函数值
return h * sum # 乘以步长得到近似积分值
# 使用复合梯形公式计算积分
a = 1
b = 2
n = 10
trapezoidal_integral = composite_trapezoidal_rule(a, b, n)
print("复合梯形公式计算的积分值:", trapezoidal_integral)
&&&
def composite_simpson_rule(a, b, n):
# a 和 b 是积分区间的端点，n 是小区间的数量，要求 n 为奇数
h = (b - a) / n # 计算步长
sum = f(a) + f(b) # 端点的函数值
for i in range(1, n - 1):
if i % 2 == 0: # 偶数区间使用辛普森公式
sum += 2 * f(a + i * h)
else: # 奇数区间使用梯形公式
sum += 4 * f(a + i * h)
return h / 3 * sum # 乘以步长和系数1/3得到近似积分值
# 使用复合辛普森公式计算积分
n = 9 # 确保n为奇数
simpson_integral = composite_simpson_rule(a, b, n)
print("复合辛普森公式计算的积分值:", simpson_integral)
&&&
import numpy as np
import matplotlib.pyplot as plt
# 定义被积函数
def f(x):
return np.exp(-x)
# 计算三点高斯求积的权重和高斯点
def get_gauss_points_weights(n=3):
# 对于n=3的高斯求积，高斯点和权重是固定的
if n == 3:
x = [0, np.sqrt(3)/3, 2*np.sqrt(3)/3]
w = [5/9, 8/9, 5/9]
else:
raise ValueError("Only n=3 is supported for this example.")
return x, w
# 三点高斯求积公式
def gauss_quadrature_integral(a, b, n, x, w):
integral = 0
for i in range(n):
integral += w[i] * f((a * (1 - x[i]) + b * (x[i])))
return 0.5 * (b - a) * integral # 由于高斯点是在[0, 1]区间的，需要乘以区间宽度
# 将区间20等分
n_subintervals = 20
subinterval_width = 1 / n_subintervals
# 高斯点和权重
x, w = get_gauss_points_weights()
# 计算每个子区间的积分，然后求和
total_integral = 0
for i in range(n_subintervals):
# 计算当前子区间的端点
a_subinterval = i * subinterval_width
b_subinterval = (i + 1) * subinterval_width
# 计算当前子区间的积分
total_integral += gauss_quadrature_integral(a_subinterval, b_subinterval, 3, x, w)
print("使用三点高斯复合求积公式计算的积分值:", total_integral)
# 为了可视化，我们可以绘制被积函数和近似的积分区域
x_values = np.linspace(0, 1, 400)
y_values = f(x_values)
plt.plot(x_values, y_values, label='Function f(x)')
plt.scatter([a_subinterval, b_subinterval], [f(a_subinterval), f(b_subinterval)], color='red', zorder=5)
plt.fill_between(x_values, y_values, alpha=0.3)
plt.legend()
plt.show()
扫一扫
评论
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
查看更多评论
添加红包
