隐方程求解一维抛物型方程（热传导方程）

记录一下数值计算课大作业

设时间t和空间x

离散：x${}_i$=ih , t${}_n$=n$\Delta$*t (i = 0,1,…I , n = 0,1,…)
设 U（x${}_i$，t${}_n$），在该点：
对t用向前差商：$\frac{\partial u}{\partial t}$ $\approx$ $\frac{U_i^{n+1}-U_i^n}{\Delta t}$
对用中心二阶差商：$\frac{{\partial }^{2}u}{\partial x^2}$ $\approx$ $\frac{U_{i+1}^n-2U_i^n+U_{i-1}^n}{h^2}$

代入上式，得到隐格式：
$\frac{U_i^{n+1}-U_i^n}{\Delta t}$ – b*$\frac{U_{i+1}^n-2U_i^n+U_{i-1}^n}{h^2}$ = 0

化简得：–$\frac{b\ast \Delta t}{h^2}\ast U_{i-1}^{n+1}+(1+\frac{2b\ast \Delta t}{h^2})\ast U_i^{(n+1)}-\frac{b\ast \Delta t}{h^2}\ast U_{i-1}^n=U_i^n$

写出当n=0.时i=0,1,…I的式子，利用初值和边值条件，可得系数矩阵：

$$A=\left[\begin{array}{cccccccc}1+\frac{2b\ast \Delta t}{h^2}& -\frac{b\ast \Delta t}{h^2}& & & & & & \\ -\frac{b\ast \Delta t}{h^2}& 1+\frac{2b\ast \Delta t}{h^2}& -\frac{b\ast \Delta t}{h^2}& & & & & \\ & -\frac{b\ast \Delta t}{h^2}& 1+\frac{2b\ast \Delta t}{h^2}& & & & & \\ & ...& ...& ...& & & & \\ & & ...& ...& ...& & & \\ & & & & & & & \\ & & & ...& ...& ...& & \\ & & & & ...& ...& ...& \\ & & & & & -\frac{b\ast \Delta t}{h^2}& 1+\frac{2b\ast \Delta t}{h^2}& \end{array}\right]$$

令AU=C,其中：
U$=(U_1^0,U_2^0,...U_{I-1}^1)$
C$=(\frac{b\ast \Delta t}{h^2}\ast U_0^1+U_1^0,U_2^0,...U_{I-2}^0+\frac{b\ast \Delta t}{h^2}\ast U_{I-1}^0)$

可解出U，同理只需解出当n=1，2，…，时的U即可。
用matlab实现：

function sol = HeatExp2(h,maxt,deltat,a)  %a为方程中的系数b，可任取
I = 1/h; 
N = maxt/deltat;
sol = zeros(fix(I+1),fix(N+1));       
mu = a*deltat/h^2;
vetx = [0:h:1];       
veti = ones(I-1,1)    
vetii = ones(I-2,1)   
A=sparse(diag((1+2*mu)*veti,0)-diag(vetii*mu,1)-diag(vetii*mu,-1))
[L U]=lu(A)
for i=2:ceil((I+1)/2) 
    sol(i,1)=2*vetx(i);
    sol(I+2-i,1)=sol(i,1);
end
b = sol(2:I,1)        
for j=2:N             
    u=U\(L\b);        
    for i=2:I   
        sol(i,j)= u(i-1);      
    end
    b = u;
    %分别画出初始、中间和最末时刻的图：            
    if j==2           
        plot([0:h:1],sol(:,j),'c')
        axis([0 1 0 1]);
        xlabel('x');
        ylabel('u');
        hold on;
    elseif j==N/2
        plot([0:h:1],sol(:,j),'r')
        axis([0 1 0 1]);
        xlabel('x');
        ylabel('u');
        hold on;   
    elseif j==N
        plot([0:h:1],sol(:,j),'y')
        axis([0 1 0 1]);
        xlabel('x');
        ylabel('u');
        hold on;
    end
    legend('t=0.001s','t=0.05s','t=0.1s');
end

在命令行输入：

得到结果矩阵：

得到二维图像：

画出三维图像：

figure(2)
mesh([0:h:1],[0:deltat:maxt],sol)

计算结果是合理的，因为热量逐渐扩散，至端点处消失。