东南大学数值分析第七章上机作业Python实现

正文开始，题目如下：

程序设计思路如下，采用7.1.37的矩阵形式进行求解，两个系数方阵都是严格对角占优的，所以一定可逆，为了编程简单没有采用迭代格式，直接使用矩阵求逆法求解方程组。

下图为第二问程序；第三问略，具体思路如下：步长二分即M或N增倍。

import numpy as np
from sympy import *

# 此函数为调参入口
def u(x,t):
    if t==0:
        u = math.exp(x)
    elif x==0:
        u = math.exp(t)
    elif x==1:
        u = math.exp(1+t)
    return u

# 此函数为调参入口
def f(x,t):
    f = 0
    return f

a = 1
M = 640         # 此处为调参入口
N = 640         # 此处为调参入口
h = 1/M        # 此处默认x取0-1
tao = 1/N      # 此处默认t取0-1，也即T=1
x_i = np.array([i for i in np.arange(0, (M+5)*h, h)])  # (M+4)*1,包含0,索引为0~M+4,索引0和M为边界(x=0和x=1)
t_i = np.array([i for i in np.arange(0, (N+5)*tao, tao)])  # (M+4)*1,包含0,索引为0~M+4,索引0和M为边界（t=0和t=1）
r = a*tao/(h*h)
A_1 = np.zeros([M-1,M-1])    # （M-1）*（M-1）
A_2 = np.zeros([M-1,M-1])
for i in range(1,M-2):
    A_1[0,0]=1+r
    A_1[0,1]=-r/2
    A_1[M-2,M-2]=1+r
    A_1[M-2,M-3]=-r/2
    A_1[i,i]=1+r
    A_1[i,i+1]= -r/2
    A_1[i,i-1]= -r/2

    A_2[0, 0] = 1 - r
    A_2[0, 1] = r / 2
    A_2[M - 2, M - 2] = 1 - r
    A_2[M - 2, M - 3] = r / 2
    A_2[i,i]=1-r
    A_2[i, i - 1] = r/2
    A_2[i, i + 1] = r/2
uik = np.zeros([M+1,M-1])  # 共41行（t=0-1，间隔为0.025）,39列，每行t不变；
uik[0,:] = [u(x,0) for x in np.arange(h,1,h)]  # uik为解，x属于h~39*h，第i行为t = i*tao解
for i in range(1,uik.shape[0]):      # i=1,2,...,39,40(对应t=1）
    U = np.zeros(M-1)      # (M-1)*1
    U[0] = tao*f(x_i[1],t_i[i]+tao/2) + r/2*(u(0,t_i[i])+u(0,t_i[i+1]))
    U[M-2] = tao*f(x_i[M-1],t_i[i]+tao/2) + r/2*(u(1,t_i[i])+u(1,t_i[i+1]))
    for j in range(1,M-2):
        U[j] = tao*f(x_i[j+1],t_i[i]+tao/2)
    uik[i,:] = np.dot(np.dot(np.linalg.inv(A_1),A_2),uik[i-1,:])+np.dot(np.linalg.inv(A_1),U)
# print(uik)
print("u(0.2, 1.0) = ", uik[int(1/tao),int(0.2/h-1)])
print("u(0.4, 1.0) = ", uik[int(1/tao),int(0.4/h-1)])
print("u(0.6, 1.0) = ", uik[int(1/tao),int(0.6/h-1)])
print("u(0.8, 1.0) = ", uik[int(1/tao),int(0.8/h-1)])