【Matlab】抛物方程Crank-Nicolson格式

1 题目

(1) 编制用Crank-Nicolson格式求抛物方程数值解的通用程序。
$$\{\begin{array}{ll}\frac{\partial u}{\partial t}-a\frac{{\partial }^{2}u}{\partial x^2}=f(x,t)& (0<x<1,0\le t\le \mathrm{T})\\ u(x,0)=\phi (x)& (0\le x\le 1)\\ u(0,t)=\alpha (t),u(1,t)=\beta (t)& (0<t\le T)\end{array}$$

(2) 就$a=1$,$f(x,t)=0$,$\varphi (x)=exp(x)$,$\alpha (t)=exp(t)$,$\beta (t)=exp(t)$,$M=40$,$N=40$,输出点$(0.2,1.0)$,$(0.4,1.0)$,$(0.6,1.0)$,$(0.8,1.0)$这四点处$u(x,t)$的数值解。

(3) 已知所给方程的精确解为$u(x,t)=exp(x+t)$，将步长反复二分，从点$(0.2,1.0)$,$(0.4,1.0)$,$(0.6,1.0)$,$(0.8,1.0)$四点处精确解与数值解的误差观察当步长缩小一半时误差以什么规律缩小。

2 程序

1 通用程序

function U = CrankNicolson(a,M,N,phi,alpha,beta)
    h = 1 / M;
    tau = 1 / N;
    r = a * tau / h / h;

    M1 = zeros(M-1, M-1);
    M2 = zeros(M-1, M-1);
    M3 = zeros(M-1, 1);
    U = zeros(N+1, M-1);

    for i=1:M-1
        M1(i,i) = 1 + r;
        M2(i,i) = 1 - r;
        U(1,i) = phi(tau*i);
        if i-1 >= 1
            M1(i, i-1) = - r / 2;
            M2(i, i-1) = r / 2;
        end
        if i+1 <= M-1
            M1(i, i+1) = -r / 2;
            M2(i, i+1) = r / 2;
        end
    end

    for k = 1:N
        M3(1) = r /2 * ( alpha(tau*(k-1)) + alpha(tau*k) );
        M3(end) = r /2 * ( beta(tau*(k-1)) + beta(tau*k) );
        temp = inv(M1) *  M2 * U(k,:)' + inv(M1) * M3 ;
        U(k+1,:) = temp';
    end
end

2 题目二主函数

clc;clear all

%% 初始数据
a = 1;
M = 40;
N = 40; 

syms x t
phi(x) = exp(x);
alpha(t) = exp(t);
beta(t) = exp(1+t);

U = CrankNicolson(a,M,N,phi,alpha,beta);

for j = 0.2:0.2:0.8
    temp = U(end, round(j * M));
    fprintf("输出点（%0.1f,1.0)的数值解为：%f\n",j, temp)
end

3 题目三主函数

clc;clear all

%% 初始数据
a = 1;

syms x t
phi(x) = exp(x);
alpha(t) = exp(t);
beta(t) = exp(1+t);

k = 5;
step = 5;

y = zeros(1,4);
e = zeros(k,4);
for j = 1:4
    y(1,j) = exp(0.2 * j + 1);
end

for i=1:k
    M = step * 2 ^(i-1);
    N = step * 2 ^(i-1); 
    U = CrankNicolson(a,M,N,phi,alpha,beta);
    fprintf("当空间步长为%0.3f时：\n",1 / M)
    for j = 1:4
        temp = U(end, round(j * M / 5));
        e(i,j) = abs(temp - y(1,j));
        fprintf("输出点（%0.1f,1.0) = %f，误差为：%0.8f\n",j, temp,e(i,j))
    end
    fprintf("\n")
    
end

3 运行结果

题目2结果

题目3结果