隐式差分+追赶法求解PDE

本文主要是介绍隐式差分+追赶法求解PDE，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

前言

偏微分方程的求解一般分为两步：

（1）利用有限差分构造出三对角矩阵，而有限差分又分为三种：显式差分、隐式差分、C-N差分（即六点差分）有限差分学习笔记-CSDN博客

（2）追赶法求解三对角矩阵：

         三对角矩阵算法（英语：tridiagonal matrix algorithm），又称为托马斯算法（Thomas algorithm，名称源于英国数学家卢埃林·托马斯）是数值线性代数中的一种算法，通过简化形式的高斯消元法求解三对角矩阵。包含n个未知数的三对角方程组可以写成

其中a1=0、 cn=0写成矩阵形式则为

高斯消元法在求解一般线性方程组时需要时间复杂度，但对于三对角系统则只需复杂度。

原理

三对角矩阵算法可分为如下两步进行。第一步求解系数和：

以及

第二步通过回代得到最终结果：

举例

实现

clear;%清除工作区变量
clc;%清屏
close all;%关闭所有图形窗口z=[];
h=8.6125; %空气交换系数
%% 材料参数输入
m1=6;m2=60;m3=36;m4=50;% 分别对四种介质分割
m=m1+m2+m3+m4;% 介质分割和
n=5400;% 对时间分割
t=5400;% 总时长
l1=0.6/1000;l2=6/1000;l3=3.6/1000;l4=5/1000;% 四种材料厚度
lam_1=0.082;lam_2=0.37;lam_3=0.045;lam_4=0.028;% 四种材料的热传导率
de_1=300;de_2=862;de_3=74.2;de_4=1.18;% 四种材料的密度
c1=1377;c2=2100;c3=1726;c4=1005;% 四种材料的比热容%% 计算热扩散率
a1=lam_1/(c1*de_1);% I层材料的热扩散率
a2=lam_2/(c2*de_2);% II层材料的热扩散率
a3=lam_3/(c3*de_3);% III层材料的热扩散率
a4=lam_4/(c4*de_4);% IV层材料的热扩散率%% 材料长度分割和时间步长分割求解
derta_x1=l1/m1;% I层材料的分割长度
derta_x2=l2/m2;% II层材料的分割长度
derta_x3=l3/m3;% III层材料的分割长度
derta_x4=l4/m4;% IV层材料的分割长度
derta_t=t/n;% 时间步长分割%% 计算各层介质剖分的步长比
r1=derta_t/derta_x1^2*a1;% 第I层介质剖分的步长比
r2=derta_t/derta_x2^2*a2;% 第II层介质剖分的步长比
r3=derta_t/derta_x3^2*a3;% 第III层介质剖分的步长比
r4=derta_t/derta_x4^2*a4;% 第IV层介质剖分的步长比u=zeros(m+1,n+1);% 定义四层耦合介质温度分布矩阵
%% 初始条件和边界条件
u(:,1)=37;%初始条件
u(1,:)=75;%边界条件%% 差分格式的系数矩阵的构造
A=zeros(m,m);
for i=1:m1-1A(i,i)=1+2*r1;A(i,i+1)=-r1;if i>=2A(i,i-1)=-r1;end
end
A(m1,m1)=(lam_1/derta_x1+lam_2/derta_x2);
A(m1,m1-1)=-lam_1/derta_x1;
A(m1,m1+1)=-lam_2/derta_x2;for i=m1+1:m1+m2-1A(i,i)=1+2*r2;A(i,i+1)=-r2;A(i,i-1)=-r2;
end
A(m1+m2,m1+m2)=(lam_2/derta_x2+lam_3/derta_x3);
A(m1+m2,m1+m2-1)=-lam_2/derta_x2;
A(m1+m2,m1+m2+1)=-lam_3/derta_x3;for i=m1+m2+1:m1+m2+m3-1A(i,i)=1+2*r3;A(i,i+1)=-r3;A(i,i-1)=-r3;
end
A(m1+m2+m3,m1+m2+m3)=(lam_3/derta_x3+lam_4/derta_x4);
A(m1+m2+m3,m1+m2+m3-1)=-lam_3/derta_x3;
A(m1+m2+m3,m1+m2+m3+1)=-lam_4/derta_x4;for i=m1+m2+m3+1:m1+m2+m3+m4-1A(i,i)=1+2*r4;A(i,i-1)=-r4;A(i,i+1)=-r4;
end
A(m,m)=h+lam_4/derta_x4;
A(m,m-1)=-lam_4/derta_x4;%% 构造右端项
for k=2:n+1d=zeros(m,1);for i=2:m-1d(i,1)=u(i+1,k-1);endd(1,1)=u(2,k-1)+r1*u(1,k);d(m1,1)=0;d(m1+m2,1)=0;d(m1+m2+m3,1)=0;d(m,1)=37*h;%% 追赶法求解b=diag(A);a=[0;diag(A,-1)];c=diag(A,1);N=length(b);L=zeros(N);L(1)=b(1);D(1)=d(1)/L(1);C(1)=c(1)/L(1);for i=2:(N-1)L(i)=b(i)-a(i)*C(i-1);D(i)=(d(i)-D(i-1)*a(i))/L(i);C(i)=c(i)/L(i);endL(N)=b(N)-a(N)*C(N-1);D(N)=(d(N)-D(N-1)*a(N))/L(N);x(N)=D(N);for i=(N-1):-1:1x(i)=D(i)-C(i)*x(i+1);endu(2:m+1,k)=x';
end%% 绘制不同时刻不同厚度温度分布图
x=1:1:m+1;
t=1:1:t+1;
surf(t,x,u)
shading interp
xlabel('t')
ylabel('x')
zlabel('u')u=round(u,2);
xlswrite('不同时间不同厚度下的温度分布.xlsx',u)% 生成温度分布的excle文件%% 四个临界面下的部分温度分布表
U=zeros(5401,4);
U(:,1)=u(m1+1,:)';
U(:,2)=u(m1+m2+1,:)';
U(:,3)=u(m1+m2+m3+1,:)';
U(:,4)=u(m1+m2+m3+m4+1,:)';xlswrite('problem1.xlsx',U)% 存储生成四个临界面下的部分温度分布表于problem1figure
subplot(2,2,1)
plot(U(:,1),'r')
xlabel('t');ylabel('u');
title('临界面I')
axis([0 5400 30 80])
subplot(2,2,2)
plot(U(:,2),'r')
xlabel('t');ylabel('u');
title('临界面II')
axis([0 5400 30 80])
subplot(2,2,3)
plot(U(:,3),'r')
title('临界面III')
axis([0 5400 30 80])
xlabel('t');ylabel('u');
subplot(2,2,4)
plot(U(:,4),'r')
title('临界面IV')
axis([0 5400 30 80])
xlabel('t');ylabel('u');

这篇关于隐式差分+追赶法求解PDE的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！

---