追赶法/托马斯算法在常微分方程两点边值问题中的差分求解、紧差分求解（用 python 代码实现并作图）

偏微分数值解法课程作业系列

Chapter 1 常微分方程两点边值问题的差分解法

题目引入：

也可为
注：差分法紧差分法最终实现形式都为追赶法求解 代码通用

分析：

6题分析

第7题紧差分法分析

下附上第七题紧差分法代码实现：

import matplotlib.pyplot as plt
import numpy as np

def thomas(N,a,b,c,d):
    l = np.zeros((N+1,1))
    u = np.zeros((N+1,1))
    x = np.zeros((N+1,1))
    y = np.zeros((N+1,1))
    u[0] = b[0]
    for i in range(N):
        l[i+1] = a[i+1]/u[i]
        u[i+1] = b[i+1] - c[i]*l[i+1]
        y[i+1] = d[i+1] - l[i+1]*y[i]
    
    y[0] = d[0]

    x[N] = y[N]/u[N]
    
    for i in range(N):
        x[N-i-1] = (d[N-i-1]-c[N-i-1]*x[N-i])/u[N-i-1]
    
    return x

def f_x(x:float):
    return ((x**4)/20 - 6)*x

def exact_solve(x):
    return x**5/20+x**3

if __name__ == "__main__":
    #print("please input the value of N:")
    N = int(4)
    #print("thx\n")
    h = 1/N

    a = [0.0,-1+h**2/12,-1+h**2/12,-1+h**2/12,0.0]

    b = [1.0,2+10*h**2/12,2+10*h**2/12,2+10*h**2/12,1.0]

    c = [0.0,-1+h**2/12,-1+h**2/12,-1+h**2/12]

    d = []
    for i in range(N+1):
        if (i==0):
            d.append(0.0)
        elif(i==N):
            d.append(21/20)
        else:
            value = f_x((i-1)*h)+10*f_x((i)*h)+f_x((i+1)*h)
            d.append(h**2*value/12)


    print("a=",a)
    print("b=",b)
    print("c=",c)
    print("d=",d)
    u = thomas(N,a,b,c,d)
    print("\n u=",u)
    exact_solution = []
    for i in range(N+1):
        exact_solution.append(exact_solve(i*h))

    print("\n exact_solve:",exact_solution)
    
    x = [0,h,2*h,3*h,4*h]
    y_u = [0,0.06158,0.22156,0.53750,1.05]
    y_exact = [0,0.01567,0.12656,0.43374,1.05]
    plt.plot(x,y_u,color='blue',linestyle='dashed',marker='o',label="caculate solution")
    plt.plot(x,y_exact,color='red',linestyle='dashed',marker='o',label="exact solution")
    plt.legend()
    plt.show()

Reference:

[1] blog1
[2] blog2