cn.m

C-N解偏微分方程的程序。解的是定步长的抛物型偏微分方程： du/dx - a * d2u/dx2 = 0 在程序中可以更改 a 的值以实现不同系数的解。在循环中改变 a 的值以实现变系数

% Result 是数值解，h 是空间步长，t是时间步长，a是系数
clear all
format long ;    % 设置输出数据的格式

h=1/10 ; t=1/10; a=1;

r=a*t/h^2;         % 步长比

% 0<=x<=1时的初值条件的数值解(共1/h+1个点)
i=1:1/h+1;
x1=0+(i-1)*h;
Result(1,i)=exp(x1);

% 0<t<=1时的边界条件的数值解(共1/t个点)
j=2:1/t+1 ;
x2(j-1)=0+(j-1)*t;
Result(j,1)=exp(x2(j-1)); 
Result(j,1/h+1)=exp(1+x2(j-1));

e = ones(1/h-1,1);
A = spdiags([e*(-r/2) e*(1+r) e*(-r/2)], -1:1, 1/h-1, 1/h-1);   % A 矩阵的生成
B = spdiags([e*(r/2) e*(1-r) e*(r/2)], -1:1, 1/h-1, 1/h-1);   % B 矩阵的生成

% 开始计算其他结点的数值解
n=size(A,1);
for TimeIndex=1:n+1
    
    fun(1,1)=r/2*(Result(TimeIndex,1)+Result(TimeIndex+1,1));
    for vv=2:1/h-2         % 函数项
        fun(vv,1)=0;
    end
    fun(1/h-1,1)=r/2*(Result(TimeIndex,1/h+1)+Result(TimeIndex+1,1/h+1));
    
    %% 计算方程右边常数项的数值向量
    u1=Result(TimeIndex,2:n+1);
    b=B*u1'+fun;  % 三对角方程右边数据
    Result(TimeIndex+1,2:1/t) = A\b ;
end

Result ;