一维热传导方程的求解

隐式后向欧拉法求解一维抛物线方程 18377*** WingFlyM

1.求解原理描述

$<\begin{bmatrix}1+2r&-r&0& &0&0\\-r&1+2r&-r& &0&0\\0&-r&1+2r&{\cdots}&0&0\\{\vdots}&{\vdots}&{\vdots}&{\ddots}\\{0}&{0}&0&{\cdots}&{1+2r}&-r\\{0}&{0}&0&{\cdots}&0&{1+2r}\\\end{bmatrix}\begin{bmatrix}u_1^k\\u_2^k\\u_3^k\\\vdots\\u_{M-2}^k\\u_{M-1}^k\end{bmatrix}=\begin{bmatrix}u_1^k+ru_0^k\\u_2^{k-1}\\u_3^{k-1}\\\vdots\\u_{M-2}^{k-1}\\u_{M-1}^{k-1}+ru_{M}^k\end{bmatrix}>$

2.求解程序的MATLAB实现

%隐式后向欧拉法求解一维抛物线方程的MATLAB实现

function [u,x,t] = Yinshi(A,xf,T,it0,bx0,bxf,M,N)
%求解方程 u_xx = u_t ,0 <= x <=xf,0 <= t <= T
%初值： u(x,0) = it0(x)
%边界条件： u(0,t) = bx0(t)
%xf为x的取值下限
%it0 为在边界 t = 0上的函数值
%bx0为在边界x = 0上的函数值
%bxf为在边界 x= xf上的函数值
% M ---沿x轴的等分段数
% N ---沿t轴的等分段数

dx = xf/M; x=[0:M]'*dx;
dt = T/N; t = [0:N]*dt:

for i=1:M+1
   u(i,1) = it0(x(i)); 
end

for n=1:N+1
    u([1 M+1],n) = [bx0(t(n)); bxf(t(n))] 
end

r = A*dt/dx/dx;
r2 =1+2*r;
for i =1:M-1
    P(i,i) = r2;
    if i>1
        P(i-1,i) = -r;
        P(i,i-1) = -r;
    end
end


for k = 2:N+1
    b = [r*u(1,k);zeros(M-3,1); r*u(M+1,k)];
    u(2:M,k) = linsolve(P,b);
end

%设定初值求解热传导方程（典型的一维抛物线偏微分方程）

function heat_solve()
A =0.5; %热传导方程的系数
it0 = inline('sin(pi*x)','x'); %初始条件
bx0 = inline('0');
bxf = inline('0'); %边界条件
xf = 2;
M = 80;
T = 0.1;
N = 100;
[u1,x,t] = Yinshi(A,xf,T,it0,bx0,bxf,M,N);
mesh(t,x,u1); grid on
xlabel('x')
ylabel('y')
zlabel('U')