二维热导方程Matlab数值解案例

过冷水之前有和大家分享一维热传导方程的有限差分法求解的具体实现过程。本次就和大家一起来看看如何根据一维热传导有限差分法的思想求解二维热传导方程进行求解。

形式如下的二维热传导方程：

边界条件是：

初值为：

差微分方法思路：

变形为迭代式如下：

    根据该迭代过程就可求出任意时刻的温度T随空间的分布情况。现在回到我们之前和大家分享的二维微分方程的具体案例中看一下怎么解：

clear all
a=0.5;
u_xy0=inline('0','x','y');
u_xyt=inline('x^2*cos(pi*y)-y^2*cos(pi*x)','x','y','t');
f=inline('0','x','y','t');
D=[0,2*pi,0,2*pi];
T=500;
Mx=100;
My=100;
N=100;
ox=(D(2)-D(1))/Mx;
x=D(1)+[0:Mx]*ox;
oy=(D(4)-D(3))/My;
y=D(3)+[0:My]*oy;
ot=T/N;
t=[0:N]*ot;
%初始化 u
for j=1:Mx+1;
    for i=1:My+1;
        u(i,j)=u_xy0(x(j),y(i));
    end
end
rx=a*ot/(ox*ox);
ry=a*ot/(oy*oy);
rx1=1+2*rx;
rx2=1-2*rx;
ry1=1+2*ry;
ry2=1-2*ry;
for j=1:Mx-1;
    A(j,j)=ry1;
    if j>1
        A(j-1,j)=-ry;
        A(j,j-1)=-ry;
    end
end
%A为y方向隐式时的系数矩阵
for i=1:My-1;
    B(i,i)=rx1;
    if i>1;
        B(i-1,i)=-rx;
        B(i,i-1)=-rx;
    end
end
%B为x方向隐式时的系数矩阵
for k=1:N
    t=k*ot;u_1=u;
    %for m=1:Mx+1;


        %for n=1:My+1;
            %v(n,m)=feval(f,x(m),y(n),t);
        %end
    %end
    u_1=u;
    for i=1:My+1
        u(i,1)=feval(u_xyt,x(1),y(i),t);
        u(i,Mx+1)=feval(u_xyt,x(Mx+1),y(i),t);
    end
    for j=1:My+1;
        u(1,j)=feval(u_xyt,x(j),y(1),t);
        u(My+1,j)=feval(u_xyt,x(j),y(My+1),t);
    end
    if mod(k,2)==0;
        for i=2:My;
            jj=2:Mx;
            bx=[ry*u(i,1),zeros(1,Mx-3),ry*u(i,My+1)]+rx*(u_1(i-1,jj)+u_1(i+1,jj))+rx2*u_1(i,jj)+ot*u_1(i,jj);
            opts.UT = true; opts.TRANSA = true;
            u(i,jj)=linsolve(A,bx',opts)';
        end
    else
        for j=2:Mx
            ii=2:My;
            by=[rx*u(1,j);zeros(My-3,1);rx*u(Mx+1,j)]+ry*(u_1(ii,j-1)+u_1(ii,j-1))+ry2*u_1(ii,j)+ot*u_1(ii,j);
            opts.UT = true; opts.TRANSA = true;
            u(ii,j)=linsolve(B,by,opts);
        end
    end
end
mesh(x,y,u);
xlabel('x');
zlabel('温度T')

过冷水本期要和大家分享的二维热传导的问题就是这些了，热传导属于比较复杂的问题，物理模型不容易理解，解决了物理模型其实就是偏微分方程的求解问题，微分方程求解的方法有很多，Matlab爱好者后期会一直分享。

往期回顾>>>>>>

热导方程的Matlab数值解方法

从泰勒级数说傅里叶级数

你所不知道的Monte Carlo形式

数值计算——MATLAB数值积分原理详讲

基于Hough变换原理实现图像直线检测【附源代码】