Matlab 偏微分方程有限差分法

计时开始


开始计时并清除所有之前变量的定义，确保结果的准确性。


```matlab




tic


clear all


```


输入参数定义

设定了关键参数包括扩散系数 `a`、结束时间节点集合 `T`、空间步长 `dx` 和时间节点在空间上的分布、时间步长 `dt` 以及计算 `λ=adt/dx/dx`。

```matlab


a = 1;


T = [0.01, 0.05, 0.1]; % 结束时间节点


dx = 0.1; % 空间步长


x = 0:dx:1; % 空间范围剖分


dt = 0.01; % 时间步长

t = 0:dt:T(end);       % 最大时间范围剖分, 充分覆盖所有时间节点

lambda = dt/dx/dx; % λ值计算


```


解析解构建

这里采用级数的方式计算了不同时间节点对应的算例二的解析解。以图的形式展示了解析解随时间的变化。

```matlab


xx = 0:0.01:1;


for tt = 1:size(T,2) % 对于三个时间节点进行解析解构造


ff = zeros(1, size(xx,2));


for i = 1:size(xx,2)


uu = 0;

for j = 1:17 % 累加17次以满足精度要求, ATTNTION: 这里的迭代次数是固定的, 基于经验或特定要求调整

f = sin(jpi/2)sin(jpixx(i))exp(jjpipiT(tt))/j^2;

uu = uu + f;


end


ff(i) = 8 uu / pi / pi;


end


subplot(1, 3, tt)


plot(xx, ff, 'b')


end


```


隐式差分格式解

利用CrankNicolson格式构造方程组解，并依随时间演进的方式计算不同时间节点的近似解。解析解与CrankNicolson格式解的对比，则通过图形直观地展示。

```matlab


n = size(x,2) 2;


A = (1 + alambda) eye(n); % 系数矩阵A


cc = 0.5alambda;


for i = 2:n


A(i,i1) = cc;


A(i1,i) = A(i,i1);


end

u = zeros(size(t,2), size(x,2)); % 初始化解矩阵

u(1,:) = 2x.(x<=0.5) + 2(1x).(x>0.5); % 初始条件运用

for j = 2:length(t) % 针对每一个时间节点进行迭代计算


B = zeros(n,1); % 初始化列向量B


for i = 2:n+1

B(i1) = 0.5alambdau(j1,i1) + (1  alambda)u(j1,i) + 0.5alambdau(j1,i+1);

end


u(j,2:n+1) = Thomas(A,B); % 使用追赶法计算未知数


end


% 矩阵赋值与绘图


for i = 1:size(T,2)


subplot(1,3,i)

plot(x,u(T(i)/dt+1,:), 'r') % 对每个时间节点的近似解绘制

str{i} = ['T=' num2str(T(i))]; % 配置图例的文本内容

str1{i} = ['S=' num2str(T(i))];  % 方程组解法名称标识

legend({str{i}, str1{i}}) % 设置并显示图例


title(['λ=', num2str(lambda)]) % 标题显示λ值


xlim([0 1])


ylim([0 1])


end


toc


```