Ⅰ 课堂练习

第1周

from scipy import linalg
import numpy as np

X_temp = [[3, 4, 5],
          [2, 5, 4],
          [4, 6, 3]]
C = [54,44,55]

# 3 * x1 + 4 * x2 + 5 * x3 = 54
# 2 * x1 + 5 * x2 + 4 * x3 = 44
# 4 * x1 + 6 * x2 + 3 * x3 = 55

X_temp = np.array(X_temp) # 系数矩阵
C = np.array(C)  # 常数列
X = linalg.solve(X_temp,C)

print(X)

[7. 2. 5.]

# 二分法
import math

def func(x):
    return math.pow(x, 3) - x - 1

a = 1
b = 1.5
count = 0

while abs(func(a)) >= 0.0001 and abs(func(b)) >= 0.0001:
    mid = (a + b) / 2
    if func(a) * func(mid) < 0:
        b = mid
    else:
        a = mid
    count += 1
    print(count,"f(a) = ",func(a),"f(b) = ",func(b))

print(count)

1 f(a) =  -0.296875 f(b) =  0.875
2 f(a) =  -0.296875 f(b) =  0.224609375
3 f(a) =  -0.051513671875 f(b) =  0.224609375
4 f(a) =  -0.051513671875 f(b) =  0.082611083984375
5 f(a) =  -0.051513671875 f(b) =  0.014575958251953125
6 f(a) =  -0.018710613250732422 f(b) =  0.014575958251953125
7 f(a) =  -0.0021279454231262207 f(b) =  0.014575958251953125
8 f(a) =  -0.0021279454231262207 f(b) =  0.006208829581737518
9 f(a) =  -0.0021279454231262207 f(b) =  0.002036650665104389
10 f(a) =  -4.659488331526518e-05 f(b) =  0.002036650665104389

# 牛顿迭代法
import math

def func(x):
    return pow(x, 3) - x - 1
def func1(x):
    return 3 * pow(x, 2) -1

a, count = 0.5, 0

while abs(func(a)) >= 0.0001:
    # x(n+1) = x(n)-f(xn) / f'(xn)
    a = a - func(a) / func1(a)
    count += 1
    print(count,"a = ", a, "f(a) = ",func(a))

1 a =  -5.0 f(a) =  -121.0
2 a =  -3.364864864864865 f(a) =  -35.73319694786094
3 a =  -2.280955053664953 f(a) =  -10.586297439073974
4 a =  -1.556276567967263 f(a) =  -3.213020231094429
5 a =  -1.0435052271790375 f(a) =  -1.0927709112857285
6 a =  -0.5614095187713111 f(a) =  -0.6155358970176101                
7 a =  -11.864344921350634 f(a) =  -1659.1926475497012
8 a =  -7.925964323903187 f(a) =  -490.9903308071819
9 a =  -5.306828631368327 f(a) =  -145.14636187274226
10 a =  -3.5682842225998948 f(a) =  -42.86543806727349
11 a =  -2.4159242097683746 f(a) =  -12.685075952386104
12 a =  -1.6476006083209067 f(a) =  -3.824955843874735
13 a =  -1.1121747148997625 f(a) =  -1.2635104429209947
14 a =  -0.6460719137732104 f(a) =  -0.6236042645542477
15 a =  1.8263234859851267 f(a) =  3.2653008379537116
16 a =  1.4637689875292796 f(a) =  0.672531206531874
17 a =  1.339865398074654 f(a) =  0.06551360110452031
18 a =  1.3249274555821327 f(a) =  0.0008936079558332644
19 a =  1.3247179981331736 f(a) =  1.743741440130009e-07

# 例：计算 115^0.5
# 令 C ^ 0.5 = x0, 即 x0^2 - C = 0
# f(xn) = xn^2 - C = 0, f'(xn) = 2xn

import math

def func(x):
    return pow(x, 2) - 115
def func1(x):
    return 2 * x

x, count = 1, 0

while abs(func(x)) >= 0.0001:
    # xn+1 = xn - f(xn)/f'(xn) = (xn + C/n)/2
    x = x - func(x) / func1(x)
    count += 1
    print(count,"xn = ", x, "f(xn) = ",func(x))

print("115^0.5 = ", x)

1 xn =  58.0 f(xn) =  3249.0
2 xn =  29.99137931034483 f(xn) =  784.4828329369799
3 xn =  16.912907246434273 f(xn) =  171.04643152648896
4 xn =  11.85622393840778 f(xn) =  25.57004607767368
5 xn =  10.777885412983993 f(xn) =  1.162813975413144
6 xn =  10.72394097347398 f(xn) =  0.002910002554045832
7 xn =  10.723805295621908 f(xn) =  1.8408471191833087e-08
115^0.5 =  10.723805295621908

# 课后作业
import math
import matplotlib.pyplot as plt
import numpy as np

def func(x):
    return np.cos(x) - 0.5

def func2(x):
    return x * x - 1

x = np.linspace(-2, 2, 1000)  # 在[-2, 2]范围内生成1000个点
y1 = func(x)
y2 = func2(x)
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 8)) # 创建画布和子图

ax1.plot(x, y1, label='cos(x) - 0.5')
ax1.axhline(0, color='black',linewidth=0.5)  # 添加横轴
ax1.axvline(0, color='black',linewidth=0.5)  # 添加纵轴
ax1.set_title('Function: cos(x) - 0.5')
ax1.set_xlabel('x')
ax1.set_ylabel('y')
ax1.legend()
ax2.plot(x, y2, label='x^2 - 1')
ax2.axhline(0, color='black',linewidth=0.5)  # 添加横轴
ax2.axvline(0, color='black',linewidth=0.5)  # 添加纵轴
ax2.set_title('Function: x^2 - 1')
ax2.set_xlabel('x')
ax2.set_ylabel('y')
ax2.legend()
plt.tight_layout() # 调整布局
plt.show()

print("\n例1：用弦截法求 cos(x) - 0.5 = 0在 [0,2] 的解")
a, b, count = 0, 2, 0

while abs(func(b)) >= 0.0001:
    # xn+1 = xn - (xn - xn-1)/[f(xn)-f(xn-1)] * f(xn)
    temp = b - (b - a) / (func(b) - func(a)) * func(b)
    a, b = b, temp
    count += 1
    print(count,"xn = ", b, "f(xn) = ",func(b))

print("\n例2：用Stewenson法求 x^2 - 1 = 0 在 0.4 附近的解")
x, count = 0.4, 0

while abs(func2(x)) >= 0.0001:
    # xn+1 = xn - f(xn)^2/[f(xn) - f(xn-f(xn))]
    x = x - func2(x)*func2(x)/(func2(x) - func2(x-func2(x)))
    count += 1
    print(count,"xn = ", x, "f(xn) = ",func2(x))

例1：用弦截法求 cos(x) - 0.5 = 0在 [0,2] 的解
1 xn =  0.7061414637186958 f(xn) =  0.26087134865098827
2 xn =  0.9929090137599724 f(xn) =  0.04625553133201865    
3 xn =  1.0547152010282026 f(xn) =  -0.006524543105986047
4 xn =  1.0470748720339365 f(xn) =  0.00010623950856891717
5 xn =  1.0471972866635268 f(xn) =  2.2909234198564832e-07

例2：用Stewenson法求 x^2 - 1 = 0 在 0.4 附近的解
1 xn =  0.9121951219512194 f(xn) =  -0.16790005948839992
2 xn =  0.9964700186380313 f(xn) =  -0.007047501955521507
3 xn =  0.9999937915702514 f(xn) =  -1.2416820952654284e-05

第2周

# 求解n元一次方程组:
# a11x1+a12x2+...+a1nxn = b1
#       a22x2+...+a2nxn = b2
#                      ...
#                 annxn = bn

import numpy as np

# 1.Gauss 顺序消去法

def gauss(A, b):
    n = len(b)
    A1 = np.column_stack([A.astype(float), b.astype(float)]) # 增广矩阵
    # 顺序消去
    for i in range(n-1):
        for j in range(i+1, n):
            A1[j] -= A1[j][i] / A1[i][i] * A1[i]
    # 回代求解
    x = np.zeros(n)
    x[n-1] = A1[n-1][n] / A1[n-1][n-1]
    for i in range(n-2, -1, -1):
        x[i] = (A1[i][n] - np.dot(A1[i][i+1:n], x[i+1:n])) / A1[i][i]
    return x

A = np.array([[3, -1, 2],
              [2, 1, -3],
              [1, 3, -2]])
b = np.array([8, -11, -3])

x = gauss(A, b)
print(x)

[0.3 1.9 4.5]

# 2.矩阵分解法

import numpy as np

# (1) LU 分解: U x = y, L y = b
# L (U x) = b, A x = b, A = L U

def LU_decomposition(A):
    n = A.shape[0]
    U = np.copy(A).astype(float)
    L = np.eye(n)
    for k in range(n-1):
        for i in range(k+1, n):
            L[i, k] = U[i, k] / U[k, k]
            U[i, k:] -= U[i, k] / U[k, k] * U[k, k:]
    return L, U

def forward_substitution(L, b):
    n = b.shape[0]
    y = np.zeros_like(b).astype(float)
    for i in range(n):
        y[i] = b[i] - np.dot(L[i, :i], y[:i])
    return y

def backward_substitution(U, y):
    n = y.shape[0]
    x = np.zeros_like(y).astype(float)
    for i in range(n-1, -1, -1):
        x[i] = (y[i] - np.dot(U[i, i+1:], x[i+1:])) / U[i, i]
    return x

def solve_lu(A, b):
    L, U = LU_decomposition(A)
    y = forward_substitution(L, b)
    x = backward_substitution(U, y)
    return x

# (2) QR 分解: 

def solve_qr(A, b):
    Q, R = np.linalg.qr(A)
    y = np.matmul(Q.T, b)
    x = np.linalg.solve(R, y)
    return x

A = np.array([[3, -1, 2],
              [2, 1, -3],
              [1, 3, -2]])
b = np.array([8, -11, -3])

x = solve_qr(A, b)
print(x)
x = solve_lu(A, b)
print(x)

[0.3 1.9 4.5]
[0.3 1.9 4.5]

# 3.迭代法

import numpy as np

# (1) 雅可比迭代法

def jacobi_iteration(A, b, x0, max_iter=100, tol=1e-6):
    n = len(b)
    x = x0.copy().astype(float)
    for k in range(max_iter):
        x_new = np.zeros_like(x)
        for i in range(n):
            x_new[i] = (b[i] - np.dot(A[i, :i], x[:i]) - np.dot(A[i, i+1:], x[i+1:])) / A[i, i]
        if np.linalg.norm(x_new - x) < tol:
            break
        x = x_new
    return x

# (2) 高斯-赛德尔迭代法

def gauss_seidel_iteration(A, b, x0, max_iter=100, tol=1e-6):
    n = len(b)
    x = x0.copy()
    for k in range(max_iter):
        for i in range(n):
            x[i] = (b[i] - np.dot(A[i, :i], x[:i]) - np.dot(A[i, i+1:], x[i+1:])) / A[i, i]
        if np.linalg.norm(A @ x - b) < tol:
            break
    return x

A = np.array([[3, -1, 2],
              [2, 1, -3],
              [1, 3, -2]]).astype(float)
b = np.array([8, -11, -3]).astype(float)
x0 = np.zeros_like(b)

x = jacobi_iteration(A, b, x0)
print(x)
x = gauss_seidel_iteration(A, b, x0)
print(x)

[ 6.05765589e+32 -1.07775160e+34 -7.34990508e+33]
[ 1.99045941e+74 -1.76836022e+75 -2.55301736e+75]

第3周

# 向前差分法
from math import *

# 计算 f(x)=x^2+e^x+log(x)+sin(x) 在 x=0.5 处的微分
def f(x):
    return x**2 + exp(x) + log(x) + sin(x)

# 步长取10^(-2)~10^(-14)
for i in range(2, 15):
    h = 10**(-i)
    df = (f(0.5 + h) - f(0.5)) / h
    print("{:.0e}, {:.10f}".format(h, df))

1e-02, 5.5224259821
1e-03, 5.5258912717
1e-04, 5.5262623253
1e-05, 5.5262996793
1e-06, 5.5263034173
1e-07, 5.5263037901
1e-08, 5.5263038812
1e-09, 5.5263038590
1e-10, 5.5263038590
1e-11, 5.5263127408
1e-12, 5.5262461274
1e-13, 5.5311311087
1e-14, 5.5511151231

# 用 Gauss 法计算 sqrt(x+1.5) 在 [-1, 1] 的积分
import numpy as np

def f(x):
    return np.sqrt(x + 1.5)

def gauss(n):
    a, b = -1, 1
    # 计算 Legendre 多项式的积分点和权重
    x, w = np.polynomial.legendre.leggauss(n)
    # 计算积分点的值
    x = 0.5 * (b - a) * x + 0.5 * (b + a)
    # 计算权重函数下的积分值
    return np.sum(w * f(x)) * 0.5 * (b - a)

print(gauss(10))

2.399529123117916

第4周

1、设降落伞从跳伞塔落下，所受空气阻力与速度成正比，初始时刻速度为0，求下落速度与时间的函数关系.

mg - kv = m dv/dt
dt/m = dv/(mg - kv)
t/m = -ln(mg - kv)/k
-kt/m = ln(mg - kv) + C
代入初始条件v(0) = 0得：C = -ln(mg)
-kt/m = ln(1 - kv/mg)
exp(-kt/m) = 1 - kv/mg
mg exp(-kt/m) = mg - kv
kv = mg(1 - exp(-kt/m))
v = (mg/k)(1 - exp(-kt/m))


2、求以下方程的通解：y' = (y+xlnx)/x
y' - y/x = lnx, P = -1/x, Q = lnx
y = exp(∫1/x)(C+∫lnx exp(-∫1/x) dx)
y = x(C+∫(lnx)/x dx)
y = x(C+∫(lnx) d(lnx))
y = x(C+(lnx)^2/2)

# Mathematica:

# 我们从一个简单的例子开始，即求解df/dx = af(x)。有两种方法可以让“DSolve”返回解，这里我们将遵循这两种方法。
# 通过练习，并根据情况，你可能更喜欢其中一个。第一种解决方法是输入。
# DSolve [f'[x] = af [x], f, x]
# 求解微分方程
# {{f→ Function[{x}, e^ax c1]}}
# 第二种解决方案在第二个参数中使用“f[x]”，即:
# DSolve[f'[x] == af[x],f[x], x]
# 求解微分方程
# {{f[x]→ exp(ax) c1}}

from sympy import Function, dsolve, Derivative, Eq, symbols

# 定义符号变量和函数
x = symbols('x')
a = symbols('a')
f = Function('f')(x)

# 定义微分方程
diff_eq = Eq(Derivative(f, x), a*f)

# 求解微分方程
print(dsolve(diff_eq, f))

Eq(f(x), C1*exp(a*x))

# Mathematica:
# 使用“f[x]”作为目标函数的“DSolve”形式返回一个更简单的，但是不太灵活的解决方案，即替换函数“f[x]”。在本书的其余部分，我将倾向于使用这种更简单的形式。
# 让我们回到替换语句这是简单微分方程的解。很明显，它用一个纯函数代替了f，这个函数就是e乘以一个常数，但是让我们试着理解这个常数的形式C[1]。
# 它是大写的，所以它一定是一个内置的Mathematica函数，Mathematica文档会告诉你“c[il]是在表示各种符号计算结果时生成的第i个参数或常数的默认形式。
# 这很清楚。你可以直接用一个值或参数替换“C[1]”。

# sol = Dsolve[f'[x] == a f[x], f[x], x];
# sol / .c[1]→b
# {{f[x]→b exp(ax)}}
# 将以下的Mathematica代码改为python代码：DSolve [f'[x] = af [x], f, x]

from sympy import Function, dsolve, Derivative, Eq, symbols

# 定义符号变量和函数
x = symbols('x')
a = symbols('a')
b = symbols('b')
f = Function('f')(x)

# 定义微分方程
diff_eq = Eq(Derivative(f, x), a*f)
# 求解微分方程
solution = dsolve(diff_eq, f)

# 设置初始条件
print(solution.subs('C1', b))

Eq(f(x), b*exp(a*x))

from sympy import Function, dsolve, Derivative, Eq, exp, symbols

x = symbols('x')
y = Function('y')(x)

# 求方程 y''-2y'-3y = 3x exp(2x) 的一个特解
diff_eq = Eq(Derivative(y, x, 2) - 2*Derivative(y, x) - 3*y, 3*x*exp(2*x))
solution = dsolve(diff_eq, y)
print(solution.rhs)

# 求方程 y''+y' = x 的一个特解
diff_eq2 = Eq(Derivative(y, x, 2) + Derivative(y, x), x)
solution2 = dsolve(diff_eq2, y)
print(solution2.rhs)

C1\*exp(-x) + C2\*exp(3\*x) + (-x - 2/3)\*exp(2\*x)
C1 + C2\*exp(-x) + x**2/2 - x

# 考虑阻力的抛体运动

# sol = Dsolve[{mv'[t] == -mg-bv[t], v(0)==v0}, v, t];
# solAlone = Part[Part[sol,1], 1];
# vel = v[t] / .solAlone;
# solt = Solve[vel = 0, t, Reals]


from sympy import symbols, Function, dsolve, solve

t, v0, mg, b = symbols('t v0 mg b')
v = Function('v')(t)
sol = dsolve(v.diff(t) + mg - b*v, v, ics={v.subs(t, 0): v0})
vel = sol.rhs
solt = solve(vel, t)
zero_points = [sol.subs(t, point) for point in solt]

print(zero_points)

[Eq(v(log(-mg/(b*v0 - mg))/b), 0)]

# 考虑阻力的抛体运动

# val = {m -> 0.1, g -> 9.8, v0 -> 20};
# Plot[t / .solt / .vals, {b, 0, 1},
# PlotRange -> {{0,0.5}, {0,2.5}}]

import matplotlib.pyplot as plt
import numpy as np

vals = {'m': 0.1, 'g': 9.8, 'v0': 20}

def calculate_t(b, m, g, v0):
    return np.sqrt(m / (b * g)) * np.arctan(np.sqrt(b * g / m) * v0)

b_values = np.linspace(0, 1, 100)
t_values = [calculate_t(b, vals['m'], vals['g'], vals['v0']) for b in b_values]

plt.plot(b_values, t_values)
plt.xlim(0, 0.5), plt.ylim(0, 2.5)
plt.xlabel('b'), plt.ylabel('t')
plt.show()

# 考虑阻力的抛体运动

# sol = Dsolve[{my'[t] == -mg-by[t], y(0)==h, y'[0]==0}, y, t];
# height = Simplify[Part[y[t]/.%, 1]]
# Limit[height, b->0]

from sympy import symbols, Function, Eq, dsolve, simplify, limit

t = symbols('t')
y = Function('y')(t)
m, g, b, h = symbols('m g b h')

diff_eq = Eq(m*y.diff(t, t), -m*g - b*y.diff(t))

# 求解微分方程
solution = dsolve(diff_eq, y)
print("Solution:", dsolve(diff_eq, y))

# 计算高度
height = solution.rhs
print("Height:",simplify(height))

# 当 b 趋近于 0 时的极限
b_limit = limit(height, b, 0)
print("Limit as b approaches 0:", b_limit)

Solution: Eq(y(t), C1 + C2\*exp(-b\*t/m) - g\*m\*t/b)
Height: C1 + C2\*exp(-b\*t/m) - g\*m\*t/b    
Limit as b approaches 0: -oo\*sign(g\*m\*t)

第5周

# 1.Euler法
# 例：用前向Euler法求解初值问题，取步长h=0.1，并与解析解 y=-1/x 比较

def forward_euler(y0, x0, xn, h):
    x_values = [x0]
    y_values = [y0]
    x = x0
    y = y0
    while x < xn:
        y += h * (-1 / x)
        x += h
        x_values.append(x)
        y_values.append(y)
    return x_values, y_values

# 初始值, 最终值, 步长
x0, y0, xn, h = 1, -1, 10, 0.1

# 使用前向 Euler 方法求解
print(forward_euler(y0, x0, xn, h))

([1, 1.1, 1.2000000000000002, 1.3000000000000003, 1.4000000000000004, 1.5000000000000004, 1.6000000000000005, 1.7000000000000006, 1.8000000000000007, 1.9000000000000008, 2.000000000000001, 2.100000000000001, 2.200000000000001, 2.300000000000001, 2.4000000000000012, 2.5000000000000013, 2.6000000000000014, 2.7000000000000015, 2.8000000000000016, 2.9000000000000017, 3.0000000000000018, 3.100000000000002, 3.200000000000002, 3.300000000000002, 3.400000000000002, 3.500000000000002, 3.6000000000000023, 3.7000000000000024, 3.8000000000000025, 3.9000000000000026, 4.000000000000003, 4.100000000000002, 4.200000000000002, 4.300000000000002, 4.400000000000001, 4.500000000000001, 4.6000000000000005, 4.7, 4.8, 4.8999999999999995, 4.999999999999999, 5.099999999999999, 5.199999999999998, 5.299999999999998, 5.399999999999998, 5.499999999999997, 5.599999999999997, 5.699999999999997, 5.799999999999996, 5.899999999999996, 5.999999999999996, 6.099999999999995, 6.199999999999995, 6.2999999999999945, 6.399999999999994, 6.499999999999994, 6.599999999999993, 6.699999999999993, 6.799999999999993, 6.899999999999992, 6.999999999999992, 7.099999999999992, 7.199999999999991, 7.299999999999991, 7.399999999999991, 7.49999999999999, 7.59999999999999, 7.6999999999999895, 7.799999999999989, 7.899999999999989, 7.9999999999999885, 8.099999999999989, 8.199999999999989, 8.299999999999988, 8.399999999999988, 8.499999999999988, 8.599999999999987, 8.699999999999987, 8.799999999999986, 8.899999999999986, 8.999999999999986, 9.099999999999985, 9.199999999999985, 9.299999999999985, 9.399999999999984, 9.499999999999984, 9.599999999999984, 9.699999999999983, 9.799999999999983, 9.899999999999983, 9.999999999999982, 10.099999999999982], [-1, -1.1, -1.190909090909091, -1.2742424242424242, -1.351165501165501, -1.4225940725940724, -1.489260739260739, -1.551760739260739, -1.6105842686725038, -1.6661398242280594, -1.7187714031754278, -1.7687714031754278, -1.8163904507944755, -1.861844996249021, -1.905323257118586, -1.9469899237852528, -1.9869899237852529, -2.0254514622467914, -2.0624884992838286, -2.098202784998114, -2.1326855436188037, -2.166018876952137, -2.198276941468266, -2.229526941468266, -2.259829971771296, -2.2892417364771784, -2.317813165048607, -2.3455909428263846, -2.372617969853412, -2.398933759327096, -2.4245747849681214, -2.4495747849681213, -2.4739650288705604, -2.497774552680084, -2.5210303666335725, -2.5437576393608454, -2.5659798615830676, -2.5877189920178503, -2.6089955877625313, -2.6298289210958647, -2.6502370843611707, -2.6702370843611707, -2.6898449274984255, -2.7090756967291947, -2.727943621257497, -2.7464621397760154, -2.7646439579578335, -2.7825011008149763, -2.800044960464099, -2.817286339774444, -2.8342354923168167, -2.8509021589834833, -2.867295601606434, -2.8834246338644984, -2.8992976497375142, -2.9149226497375142, -2.9303072651221296, -2.9454587802736447, -2.960384153407973, -2.975090035760914, -2.9895827893841025, -3.0038685036698167, -3.0179530107120702, -3.031841899600959, -3.0455405297379454, -3.059054043251459, -3.072387376584792, -3.0855452713216343, -3.0985322843086474, -3.11135279712916, -3.1240110249772615, -3.1365110249772616, -3.1488567039896074, -3.161051825940827, -3.1731000187119114, -3.1850047806166732, -3.196769486499026, -3.2083973934757704, -3.2198916463493337, -3.23125528271297, -3.2424912377691495, -3.2536023488802606, -3.2645913598692715, -3.2754609250866626, -3.2862136132587056, -3.296851911131046, -3.3073782269205196, -3.317794893587186, -3.3281041719377016, -3.3383082535703545, -3.3484092636713645, -3.3584092636713647])

# 2.Runge-Kutta法
# 例：用4阶RK法求解初值问题，取步长h=0.2, 并与解析解 y=2exp(x)-x-1 比较

import numpy as np
import matplotlib.pyplot as plt

def f(x, y):
    return 2 * np.exp(x) - x - 1

def runge_kutta_4th_order(x0, y0, h, xn):
    n = int((xn - x0) / h)
    x = np.linspace(x0, xn, n + 1)
    y = np.zeros(n + 1)
    y[0] = y0
    for i in range(n):
        k1 = h * f(x[i], y[i])
        k2 = h * f(x[i] + h/2, y[i] + k1/2)
        k3 = h * f(x[i] + h/2, y[i] + k2/2)
        k4 = h * f(x[i] + h, y[i] + k3)
        y[i+1] = y[i] + (k1 + 2*k2 + 2*k3 + k4) / 6
    return x, y

# 初始值, 最终值, 步长
x0, y0, xn, h = 0, 1, 1, 0.2

# 使用4阶Runge-Kutta方法求解
x_rk4, y_rk4 = runge_kutta_4th_order(x0, y0, h, xn)

# 解析解
x_exact = np.linspace(x0, xn, 100)
y_exact = 2 * np.exp(x_exact) - x_exact - 1

# 比较结果
plt.plot(x_rk4, y_rk4, label='RK4 Method')
plt.plot(x_exact, y_exact, label='Exact Solution')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Comparison of RK4 Method with Exact Solution')
plt.legend()
plt.grid(True)
plt.show()

# 1) Mathematica

# 函数“NDSolve”可以找到微分方程的数值解，当然，方程中不能有符号参数，
# 我们必须告诉Mathematica我们想要找到解的范围，但这很简单。对于上面的问题，我们可以尝试下面的方法:

# Remove ["Global ` *"];
# ode = D[D[f[x],x], x] == f[x]^2 - 2*x;
# ics = {f[0] == 0, f'[0] == 1};
# xmax = 10;
# sol = NDSolve[{ode,ics}, f[x], {x,0,xmax}];
# y = f[x] /. Part[Part[sol,1], 1];
# Plot[y, {x,0,xmax}]

# 2) Python

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# 定义微分方程系统
def ode_system(x, y):
    f, df = y
    return [df, f**2 - 2*x]

# 初始条件
ics = [0, 1]  # f[0] == 0, f'[0] == 1

# 解的范围
x_span = [0, 10]

# 使用solve_ivp解决ODE
solution = solve_ivp(ode_system, x_span, ics, t_eval=np.linspace(0, 10, 400))

# 解析解
x_values = solution.t
f_values = solution.y[0]

# 绘制解
plt.plot(x_values, f_values)
plt.xlabel('x')
plt.ylabel('f(x)')
plt.title('Numerical solution of the differential equation')
plt.grid(True)
plt.show()

# 1) Mathematica

# 粒子在倾斜势阱中的运动
# Remove ["Global ` *"];

# 建立方程，初始条件和最大时间：
# diffeq = D[D[f[u], u], u] == f[u]^2 - 2 f[u];
# y0a = 0.01;
# y0b = -0.99;
# y0c = -1.01;
# ica ={f'[0] == 0, f[0] == y0a);
# icb ={f'[0] == 0, f[0] == y0b);
# icc ={f'[0] == 0, f[0] == y0c);

# 解这三个初始条件的方程：
# sol = NDSolve[{diffeq, ica}, f[u], {u,0,umax}];
# ya = f[u] /. Part[Part[sol, 1], 1];
# sol = NDSolve[{diffeq, icb}, f[u], {u,0,umax}];
# yb = f[u] /. Part[Part[sol, 1], 1];
# sol = NDSolve[{diffeq, icc}, f[u], {u,0,umax}];
# yc = f[u] /. Part[Part[sol, 1], 1];

# 2) Python

import numpy as np
from scipy.integrate import odeint

# 定义微分方程
def diffeq(f, u):
    return [f[1], f[1]**2 - 2*f[0]]

# 初始条件和最大时间
y0a = [0.01, 0]
y0b = [-0.99, 0]
y0c = [-1.01, 0]
umax = 10

# 对每个初始条件求解微分方程
ua = odeint(diffeq, y0a, np.linspace(0, umax, 1000))
ya = ua[:, 0]
print("y0a 的解:", ya)

ub = odeint(diffeq, y0b, np.linspace(0, umax, 1000))
yb = ub[:, 0]
print("y0b 的解:", yb)

uc = odeint(diffeq, y0c, np.linspace(0, umax, 1000))
yc = uc[:, 0]
print("y0c 的解:", yc)

第6周

n = 25;
a = 0;
b = Pi/4;
h = (b-a) / n;
esp = 10^(-3);

y1 = y2 = z1 = z2 = Table[0,{n+1}];
x = Table[0+(i-1)*h, {i,1,n+ 1}];
t = Table[0,{1}];
t[[1]] = -0.15;

f[xxx_,yy_] := -4*yy + cos[xx];
k = 1;
y1[[-1]] = 1;

While[Abs[y1[[-1]]]] ≥ eps,
    y1[[1]] = 0;
    y2[[1]] = t[[k]];
    z1[[1]] = 0;
    z2[[1]] = 1;
    For[
        i = 1, i ≤ n,
        y1[[i+1]] = y1[[i]] + h*y2[[i]];
        y2[[i+1]] = y2[[i]] + h*N[f[x[[i]],y1[[i]]]];
        z1[[i+1]] = z1[[i]] + h*z2[[i]];
        z2[[i+1]] = z2[[i]] + h*(-4*z1[[i]]); i++];
    AppendTo[t, t[[k]] - y1[[-1]] / z1[[-1]]]; k++

ndata = Transpose[{x, y1}];

p1 = ListPlot[ndata, PlotStyle -> Red, PlotLegends -> "Numerical"];
p2 = Plot[
        1/12*(-4*Cos[2*x]+4*Cos[x]^3*Cos[2*x]-2*Sqrt[2]*Sin[2*x]+
        3*Sin[x]*Sin[2*x]+Sin[2*x]*Sin[3*x]),{x,0,Pi/4},
        PlotLegends→"Analytical"];

Show[p1, p2]

import numpy as np
import matplotlib.pyplot as plt
from math import cos, sin, sqrt

# 清除所有变量（在 Python 中不是必需的）

n = 25
a = 0
b = np.pi / 4
h = (b - a) / n
eps = 1e-3

# 初始化数组
y1 = np.zeros(n + 1)
y2 = np.zeros(n + 1)
z1 = np.zeros(n + 1)
z2 = np.zeros(n + 1)
x = np.linspace(0, np.pi / 4, n + 1)
t = [-0.15]

# 定义函数 f
def f(xx, yy):
    return -4 * yy + cos(xx)

k = 1
y1[-1] = 1

# 执行迭代
while abs(y1[-1]) >= eps:
    y1[0] = 0
    y2[0] = t[k - 1]
    z1[0] = 0
    z2[0] = 1
    
    for i in range(n):
        y1[i + 1] = y1[i] + h * y2[i]
        y2[i + 1] = y2[i] + h * f(x[i], y1[i])
        z1[i + 1] = z1[i] + h * z2[i]
        z2[i + 1] = z2[i] + h * (-4 * z1[i])
    
    t.append(t[k - 1] - y1[-1] / z1[-1])
    k += 1

# 准备绘图所需数据
ndata = np.column_stack((x, y1))

# 绘图
plt.figure()
plt.plot(ndata[:, 0], ndata[:, 1], 'r-', label='Numerical')

x_values = np.linspace(0, np.pi / 4, 100)
y_values = (1 / 12) * (-4 * np.cos(2 * x_values) + 4 * np.cos(x_values)**3 * np.cos(2 * x_values) -
                      2 * np.sqrt(2) * np.sin(2 * x_values) + 3 * np.sin(x_values) * np.sin(2 * x_values) +
                      np.sin(2 * x_values) * np.sin(3 * x_values))

plt.plot(x_values, y_values, 'b-', label='Analytical')
plt.legend()
plt.show()

第7周

pde = x*D[f[x, y], x] - 2y*D[f[x, y], y] == 0;
sol = DSolve[{pde, f[1, y] == 2y + 1}, f[x, y], {x, y}];
u = f[x, y] /. Part[Part[sol, 1], 1]
(* 1 + 2x^2 y *)

import sympy as sp

# 定义符号变量
x, y = sp.symbols('x y')
# 定义符号函数
f = sp.Function('f')(x, y)
# 定义偏微分方程
pde = x * sp.diff(f, x) - 2 * y * sp.diff(f, y)
# 求解偏微分方程
sol = sp.dsolve(pde, f)
# 定义边界条件
boundary_condition = {f.subs(x, 1): 2 * y + 1}
# 代入边界条件并提取解
u = sol.subs(boundary_condition)
# 显示解
print(u)

1 + 2\*x\*\*2*y

Ⅱ 课堂实操

实操一

(* 1. 前向欧拉 shooting，数值解，因此得出的是一个散点图 *)
n = 25;
a = 0;
b = Pi/4;
h = (b - a)/n;
eps = 10^(-3);
y1 = y2 = z1 = z2 = Table[0, {n + 1}];
x = Table[0 + (i - 1) * h, {i, 1, n + 1}];
t = {-0.15}; (* 使用列表简化赋值 *)
f[xx_, yy_] := -4 * yy + Cos[xx];
k = 1;
y1[[-1]] = 1;
While[Abs[y1[[-1]]] >= eps,
  y1[[1]] = 0;
  y2[[1]] = t[[k]];
  z1[[1]] = 0;
  z2[[1]] = 1;
  For[i = 1, i <= n, i++,
    y1[[i + 1]] = y1[[i]] + h * y2[[i]];
    y2[[i + 1]] = y2[[i]] + h * N[f[x[[i]], y1[[i]]]];
    z1[[i + 1]] = z1[[i]] + h * z2[[i]];
    z2[[i + 1]] = z2[[i]] + h * (-4 * z1[[i]]);
  ];
  AppendTo[t, t[[k]] - y1[[-1]]/z1[[-1]]];
  k++;
]
ndata = Transpose[{x, y1}];
p1 = ListPlot[ndata, PlotStyle -> Red, PlotLegends -> "Numberical"];

Show[p1]

p2 = Plot[
  1/12 * (-4 * Cos[2 * x] + 4 * Cos[x]^3 * Cos[2 * x] - 2 * Sqrt[2] * Sin[2 * x] + 3 * Sin[x] * Sin[2 * x] + Sin[2 * x] * Sin[3 * x]),
  {x, 0, Pi/4},
  PlotLegends -> "Analytical"
];
Show[p1, p2]

import numpy as np
import matplotlib.pyplot as plt

n = 25
a = 0
b = np.pi / 4
h = (b - a) / n
eps = 10**(-3)
y1 = np.zeros(n + 1)
y2 = np.zeros(n + 1)
z1 = np.zeros(n + 1)
z2 = np.zeros(n + 1)
x = np.linspace(0, b, n + 1)
t = [-0.15]  # Using a list to simplify appending
f = lambda xx, yy: -4 * yy + np.cos(xx)
k = 1
y1[-1] = 1

while np.abs(y1[-1]) >= eps:
    y1[0] = 0
    y2[0] = t[k - 1]
    z1[0] = 0
    z2[0] = 1
    for i in range(n):
        y1[i + 1] = y1[i] + h * y2[i]
        y2[i + 1] = y2[i] + h * f(x[i], y1[i])
        z1[i + 1] = z1[i] + h * z2[i]
        z2[i + 1] = z2[i] + h * (-4 * z1[i])
    t.append(t[k - 1] - y1[-1] / z1[-1])
    k += 1

ndata = np.column_stack((x, y1))
plt.plot(ndata[:, 0], ndata[:, 1], 'ro', label='Numerical')
plt.legend()

x_vals = np.linspace(0, np.pi / 4, 100)
y_vals = 1 / 12 * (-4 * np.cos(2 * x_vals) + 4 * np.cos(x_vals)**3 * np.cos(2 * x_vals) - 2 * np.sqrt(2) * np.sin(2 * x_vals) + 3 * np.sin(x_vals) * np.sin(2 * x_vals) + np.sin(2 * x_vals) * np.sin(3 * x_vals))

plt.plot(x_vals, y_vals, label='Analytical')
plt.legend()
plt.show()

[图片上传失败...(image-ab4bd-1718292605046)]

实操二

(* 2. 一阶偏微分方程解析解 *)
pde = x * D[f[x, y], x] - 2 * y * D[f[x, y], y] == 0;
sol = DSolve[{pde, f[1, y] == 2 * y + 1}, f[x, y], {x, y}];
u = f[x, y] /. Part[Part[sol, 1], 1];

import sympy as sp

x, y = sp.symbols('x y')
f = sp.Function('f')(x, y)
pde = x * sp.diff(f, x) - 2 * y * sp.diff(f, y)
sol = sp.dsolve(sp.Eq(pde, 0), f, hint='1st_exact', ics={f.subs(x, 1): 2 * y + 1})
u = sol.rhs
print(u)

1 + 2 * x ^ 2 * y

实操三

(* 3. 在数值求解微分方程组中，采用 shooting 的方法 *)
sols = Map[
  First[NDSolve[{x''[t] + Sin[x[t]] == 0, x[0] == x[10] == 0}, x, t,
    Method -> {"Shooting", "StartingInitialConditions" -> {x[0] == 0, x'[0] == #}}]] &,
  {1.5, 1.75, 2}
];
Plot[Evaluate[x[t] /. sols], {t, 0, 10}, PlotStyle -> {Black, Blue, Green}];

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# Define the ODE system
def ode_system(t, y):
    x, v = y
    dxdt = v
    dvdt = -np.sin(x)
    return [dxdt, dvdt]

# Define the shooting method function
def shooting_method(v0):
    sol = solve_ivp(ode_system, [0, 10], [0, v0], t_eval=np.linspace(0, 10, 500))
    return sol

# Solve the ODE system for different initial velocities using shooting method
initial_velocities = [1.5, 1.75, 2]
solutions = [shooting_method(v0) for v0 in initial_velocities]

# Plot the solutions
colors = ['black', 'blue', 'green']
for sol, color in zip(solutions, colors):
    plt.plot(sol.t, sol.y[0], color=color)

plt.xlabel('t')
plt.ylabel('x(t)')
plt.legend([f'v0={v0}' for v0 in initial_velocities])
plt.grid(True)
plt.show()

实操四

(* 4. 尝试使用数值解，解偏微分方程 *)
pde = D[D[u[x, y], x], x] + D[D[u[x, y], y], y] == 0;
bcs = {u[-1, y] == 0, u[1, y] == 0, u[x, -1] == 0, u[x, 1] == 10};
sol = NDSolve[{pde, bcs}, u[x, y], {x, -1, 1}, {y, -1, 1}];
z = u[x, y] /. Part[Part[sol, 1], 1];
Plot[
  {z /. x -> 0, z /. x -> 0.5, z /. x -> 0.75, z /. x -> 0.95},
  {y, -1, 1},
  PlotStyle -> {Solid, Dashed, Dotted, DotDashed},
  PlotLegends -> {0, 0.5, 0.75, 0.95}
]

import numpy as np
from scipy.integrate import solve_bvp
import matplotlib.pyplot as plt

# 定义偏微分方程
def pde(x, y, u):
    dudx = u[1]
    dudy = u[2]
    return np.vstack((dudx, dudy, np.zeros_like(x)))

# 定义边界条件
def bcs(ya, yb):
    return np.array([ya[0], ya[0], yb[0] - 10, yb[0]])

# 定义初始猜测值
x = np.linspace(-1, 1, 100)
y = np.linspace(-1, 1, 100)
u_guess = np.zeros((3, x.size))

# 解偏微分方程
sol = solve_bvp(pde, bcs, x, u_guess)

# 绘制结果图
x_plot = np.array([0, 0.5, 0.75, 0.95])
y_plot = np.linspace(-1, 1, 100)
z_plot = sol.sol(x_plot)
plt.plot(y_plot, z_plot[0], label='0', linestyle='solid')
plt.plot(y_plot, z_plot[1], label='0.5', linestyle='dashed')
plt.plot(y_plot, z_plot[2], label='0.75', linestyle='dotted')
plt.plot(y_plot, z_plot[3], label='0.95', linestyle='dashdot')
plt.xlabel('y')
plt.ylabel('u(x, y)')
plt.legend(title='x')
plt.show()

Ⅲ 课后作业

作业一

(* 1. 熟悉使用 Mathematica，求出方程的数值解，并作图 *)
amp = A;
tau = a;
osc = Sin[Pi x];
damp = Exp[-x/tau];
func1 = amp damp osc;
func2 = amp damp;
func3 = -amp damp;
Plot[{func1, func2, func3} /. {A -> 10, a -> 5}, {x, 0, 10}]

x = a Cos[q];
y = b Sin[q];
ellipse1 = {x, y} /. {a -> 3, b -> 2};
ellipse2 = {x, y} /. {a -> 3, b -> 1};
ParametricPlot[{ellipse1, ellipse2}, {q, 0, 2 Pi}]

import numpy as np
import matplotlib.pyplot as plt

A = 10
a = 5
amp = A
tau = a
x_values = np.linspace(0, 10, 1000)

osc = np.sin(np.pi * x_values)
damp = np.exp(-x_values / tau)
func1 = amp * damp * osc
func2 = amp * damp
func3 = -amp * damp

a = 3
b1 = 2
b2 = 1
q_values = np.linspace(0, 2 * np.pi, 1000)

x1 = a * np.cos(q_values)
y1 = b1 * np.sin(q_values)
x2 = a * np.cos(q_values)
y2 = b2 * np.sin(q_values)

作业二

(* 2. 熟悉偏导、积分的代码，并用其解决一些问题 *)
D[x^2, x]
D[x^a, {x, 2}]
f = Exp[-a t] Sin[Pi x];
D[f, {x, 2}, t]
Integrate[2 x, x]
Integrate[4 x y, x, y]
Integrate[Cos[x], x]
Integrate[1/x, x]
Integrate[DiracDelta[x], x]
x HeavisideTheta[x]

x = A Exp[-β t] Cos[ω t];
v = D[x, t];
a = D[v, t];
f = -k x - b v /. {k -> m (ω^2 + β^2), b -> 2 m β};
m a
a

from sympy import symbols, diff, integrate, cos, DiracDelta, Heaviside, exp, sin, pi

# Define symbols
x, a, t, y, A, beta, omega, k, m, b = symbols('x a t y A beta omega k m b')

# Differentiation
diff_x_squared = diff(x**2, x)
diff_x_to_a = diff(x**a, x, 2)
f = exp(-a * t) * sin(pi * x)
diff_f_x_2_t = diff(f, x, 2, t)

# Integration
integrate_2x = integrate(2 * x, x)
integrate_4xy = integrate(4 * x * y, (x, 0, y), (y, 0, 1))
integrate_cos_x = integrate(cos(x), x)
integrate_1_over_x = integrate(1/x, x)
integrate_dirac_delta = integrate(DiracDelta(x), x)
x_heaviside_theta = x * Heaviside(x)

# Substitution and calculations
x = A * exp(-beta * t) * cos(omega * t)
v = diff(x, t)
a = diff(v, t)
substitutions = [(k, m * (omega**2 + beta**2)), (b, 2 * m * beta)]
f_substituted = (-k * x - b * v).subs(substitutions)
m_times_a = m * a
a_result = a.subs({A: 1, beta: 0.5, omega: 1, m: 1})

Differentiation:
d(x^2)/dx: 2*x
d(x^a)/dx^2: a*x**a*(a - 1)/x**2
D[f, {x, 2}, t]: pi**2*a*exp(-a*t)*sin(pi*x)

Integration:
Integrate 2x: x**2
Integrate 4xy: 1/2
Integrate cos(x): sin(x)
Integrate 1/x: log(x)
Integrate DiracDelta(x): Heaviside(x)
x HeavisideTheta(x): x*Heaviside(x)

Substitution and calculations:
Substituted function f: -A*m*(beta**2 + omega**2)*exp(-beta*t)*cos(omega*t) - 2*beta*m*(-A*beta*exp(-beta*t)*cos(omega*t) - A*omega*exp(-beta*t)*sin(omega*t))
m * a: m*(A*beta**2*exp(-beta*t)*cos(omega*t) + 2*A*beta*omega*exp(-beta*t)*sin(omega*t) - A*omega**2*exp(-beta*t)*cos(omega*t))
Result of a: 1.0*exp(-0.5*t)*sin(t) - 0.75*exp(-0.5*t)*cos(t)

作业三

(* 3. 求解各种微分方程 *)
DSolve[f'[x] == a f[x], f[x], x]
DSolve[f[x] == a f'[x], f[x], x]
DSolve[f''[x] - 2 f'[x] - 3 f[x] == 3 x Exp[2 x], f[x], x]
DSolve[f''[x] + f'[x] == x, f[x], x]
DSolve[f[x] == a f'[x], f[x], x]

from sympy import symbols, Function, dsolve, exp

# Define symbols and function
x, a = symbols('x a')
f = Function('f')

# Differential equations
eq1 = f(x).diff(x) - a * f(x)
eq2 = f(x) - a * f(x).diff(x)
eq3 = f(x).diff(x, 2) - 2 * f(x).diff(x) - 3 * f(x) - 3 * x * exp(2 * x)
eq4 = f(x).diff(x, 2) + f(x).diff(x) - x

# Solving the differential equations
solution1 = dsolve(eq1, f(x))
solution2 = dsolve(eq2, f(x))
solution3 = dsolve(eq3, f(x))
solution4 = dsolve(eq4, f(x))

Solution 1: Eq(f(x), C1*exp(a*x))
Solution 2: Eq(f(x), C1*exp(x/a))
Solution 3: Eq(f(x), C1*exp(-x) + C2*exp(3*x) + (-x - 2/3)*exp(2*x))
Solution 4: Eq(f(x), C1 + C2*exp(-x) + x**2/2 - x)

作业四

(* 4. 求解带阻力的有垂直加速度物体的偏微分方程 *)
sol = DSolve[{m v'[t] == -m g - b v[t], v[0] == v0}, v, t];
solAlone = Part[Part[sol, 1], 1];
vel = v[t] /. solAlone;
solt = Solve[vel == 0, t, Reals]
(* 当加速度取重力加速度，其他取确切的初始值时 *)
vals = {m -> 0.1, g -> 9.8, v0 -> 20};
Plot[t /. solt /. vals, {b, 0, 1}, PlotRange -> {{0, 0.5}, {0, 2.5}}]

from sympy import symbols, Function, dsolve, solve, plot

t, v0, b = symbols('t v0 b')
v = Function('v')(t)
m, g = symbols('m g')
vals = {m: 0.1, g: 9.8, v0: 20}
diff_eq = m * v.diff(t) + m * g + b * v
sol = dsolve(diff_eq, v, ics={v.subs(t, 0): v0})
sol_alone = sol.args[1][0]
vel = sol_alone.rhs
solt = solve(vel, t)

作业五

(* 5. 求解重力加速度与时间的关系 *)
sol = DSolve[{m y''[t] == -m g - b y'[t], y[0] == h, y'[0] == 0}, y, t];
height = Simplify[Part[y[t] /. sol, 1]];
Limit[height, b -> 0]

from sympy import symbols, Function, dsolve, simplify, limit

t, h, g, b = symbols('t h g b')
y = Function('y')(t)
m = symbols('m')
diff_eq = m * y.diff(t, t) + m * g + b * y.diff(t)
sol = dsolve(diff_eq, y, ics={y.subs(t, 0): h, y.diff(t).subs(t, 0): 0})
height = simplify(sol.rhs)
limit_height = limit(height, b, 0)

(-g*m*t**2/2 + h*m)/m

作业六

(* 6. 一般的偏微分方程解 *)
ode = D[D[f[x], x], x] == f[x]^2 - 2*x;
ics = {f[0] == 0, f'[0] == 1};
xmax = 10;
sol = NDSolve[{ode, ics}, f[x], {x, 0, xmax}];
y = f[x] /. Part[Part[sol, 1], 1];
Plot[y, {x, 0, xmax}]
(* 该方程无法给出数值解，因为其不存在，只有离散点构成的图像解 *)

ode = D[f[x], x] == (1 + x)*f[x]^2 + f[x];
ics = {f[1] == -1};
xmax = 1.5;
sol = NDSolve[{ode, ics}, f[x], {x, 1, xmax}];
y = f[x] /. Part[Part[sol, 1], 1];
Plot[y, {x, 1, xmax}]

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# Define the first ODE
def ode1(x, f):
    return f[1], f[1]**2 - 2*x

# Define initial conditions and integration range for the first ODE
x_range1 = (0, 10)
f0_1 = [0, 1]
# Solve the first ODE
sol1 = solve_ivp(ode1, x_range1, f0_1, t_eval=np.linspace(*x_range1, 100))

# Define the second ODE
def ode2(x, f):
    return (1 + x) * f[1]**2 + f[1], f[0]

# Define initial conditions and integration range for the second ODE
x_range2 = (1, 1.5)
f0_2 = [-1]
# Solve the second ODE
sol2 = solve_ivp(ode2, x_range2, f0_2, t_eval=np.linspace(*x_range2, 100))

作业七

(* 7. 一个关于势阱物理图像的例子 *)
diffeq = D[D[f[u], u], u] == f[u]^2 - 2*f[u];
y0a = 0.01;
y0b = -0.99;
y0c = -1.01;
ica = {f'[0] == 0, f[0] == y0a};
icb = {f'[0] == 0, f[0] == y0b};
icc = {f'[0] == 0, f[0] == y0c};
umax = 10;
sol = NDSolve[{diffeq, ica}, f[u], {u, 0, umax}];
ya = f[u] /. Part[Part[sol, 1], 1];
sol = NDSolve[{diffeq, icb}, f[u], {u, 0, umax}];
yb = f[u] /. Part[Part[sol, 1], 1];
sol = NDSolve[{diffeq, icc}, f[u], {u, 0, umax}];
yc = f[u] /. Part[Part[sol, 1], 1];
Plot[ya, {u, 0, umax}]
Plot[yb, {u, 0, umax}]
Plot[yc, {u, 0, umax}]

(* 当考虑势阱的能量后，图像为 *)
poten = -Integrate[ξ^2 - 2*ξ, {ξ, 0, y}];
Plot[poten, {y, -1.1, 3.5}]

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# Define the differential equation
def diffeq(u, f):
    return f[1], f[1]**2 - 2*f[0]

# Define initial conditions
y0a = 0.01
y0b = -0.99
y0c = -1.01
# Define initial conditions and integration range for cases a, b, and c
ica = (0, [y0a, 0])
icb = (0, [y0b, 0])
icc = (0, [y0c, 0])
umax = 10
# Solve the differential equations for cases a, b, and c
sol_a = solve_ivp(diffeq, (0, umax), [y0a, 0], t_eval=np.linspace(0, umax, 100))
sol_b = solve_ivp(diffeq, (0, umax), [y0b, 0], t_eval=np.linspace(0, umax, 100))
sol_c = solve_ivp(diffeq, (0, umax), [y0c, 0], t_eval=np.linspace(0, umax, 100))

# Plot the potential energy function
y_vals = np.linspace(-1.1, 3.5, 100)
poten = lambda y: -np.trapz(y**2 - 2*y, y_vals)

# 例1：用弦截法计算 cos(x)-0.5=0 在区间 [0,2] 的解，要求误差小于 0.0001
import numpy as np 

def f(x): 
    return np.cos(x) - 0.5 

def secant(): 
    x0 = 0 
    x1 = 2 
    count = 0 
    while True: 
        x2 = x1 - f(x1) * (x1 - x0) / (f(x1) - f(x0)) 
        if abs(x2 - x1) < 1e-4: 
            break 
        else: 
            x0 = x1 
            x1 = x2 
            count += 1 
    print(count) 
    return x2 

ans1 = secant() 
print(ans1)

5
1.047197551205967

# 例2：用 Stewenson 方法求解 x^2-1=0 在 x0=0.4 附近的解，要求误差小于 0.0001
def stewenson(): 
    x1 = 0.4 
    count = 0 
    while True: 
        x2 = x1 - (f(x1)**2 / (f(x1) - f(x1 - f(x1)))) 
        if abs(x2 - x1) < 1e-4: 
            break 
        else: 
            x1 = x2 
            count += 1 
    print('count={}'.format(count)) 
    return x2 

ans2 = stewenson() 
print(ans2)

count=4
11.519173063162574

# 例3：牛顿迭代法 计算 115^0.5
def newton(): 
    x1 = 0.4 
    count = 0 
    while True: 
        x2 = (x1 + (115 / x1)) / 2 
        if abs(x2 - x1) < 1e-4: 
            break 
        else: 
            x1 = x2 
            count += 1 
    print('count={}'.format(count)) 
    return x2 

ans3 = newton() 
print(ans3)

count=8
10.723805294763608

# 例4：用二分法求方程 x^3-x-1＝0 在区间 [1,1.5] 内的一个实根，要求误差小于 0.0001
def bisection(): 
    x1 = 1 
    x2 = 1.5 
    count = 0 
    while True: 
        x3 = (x1 + x2) / 2 
        if f(x3) > 0: 
            x2 = x3 
        x12 = (x1 + x2) / 2 
        if f(x3) < 0: 
            x1 = x3 
        x12 = (x1 + x2) / 2 
        if count >= 12: 
            break 
        else: 
            count += 1 
    print('count={}'.format(count)) 
    return x2 

ans4 = bisection() 
print(ans4)

count=12
1.5