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

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

设时间t和空间x

离散：xi_{i}i=ih , tn_{n}n=nΔ\DeltaΔ*t (i = 0,1,…I , n = 0,1,…)
设 U（xi_{i}i，tn_{n}n），在该点：
对t用向前差商：∂u∂t\frac{\partial u}{\partial t}∂t∂u ≈\approx≈ Uin+1−UinΔt\frac{U^{n+1}_i -U^{n}_i}{\Delta t}ΔtUin+1−Uin
对用中心二阶差商：∂2u∂x2\frac{\partial^2 u}{\partial x^2}∂x2∂2u ≈\approx≈ Ui+1n−2Uin+Ui−1nh2\frac{U^{n}_{i+1} -2U^{n}_i+U^n_{i-1}}{h^2}h2Ui+1n−2Uin+Ui−1n

代入上式，得到隐格式：
Uin+1−UinΔt\frac{U^{n+1}_i -U^{n}_i}{\Delta t}ΔtUin+1−Uin – b*Ui+1n−2Uin+Ui−1nh2\frac{U^{n}_{i+1} -2U^{n}_i+U^n_{i-1}}{h^2}h2Ui+1n−2Uin+Ui−1n = 0

化简得：–b∗Δth2∗Ui−1n+1+(1+2b∗Δth2)∗Ui(n+1)−b∗Δth2∗Ui−1n=Uin\frac{b*{\Delta t}}{h^2} *U^{n+1}_{i-1}+(1+\frac{2b*\Delta t}{h^2})*U^{(n+1)}_i-\frac{b*\Delta t}{h^2}*U^n_{i-1}=U^n_ih2b∗Δt∗Ui−1n+1+(1+h22b∗Δt)∗Ui(n+1)−h2b∗Δt∗Ui−1n=Uin

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

A=[1+2b∗Δth2−b∗Δth2−b∗Δth21+2b∗Δth2−b∗Δth2−b∗Δth21+2b∗Δth2....................................−b∗Δth21+2b∗Δth2]A= \begin{bmatrix} 1+\frac{2b*\Delta t}{h^2} & -\frac{b*\Delta t}{h^2} & \\ -\frac{b*\Delta t}{h^2} & 1+\frac{2b*\Delta t}{h^2} & -\frac{b*\Delta t}{h^2} \\ & -\frac{b*\Delta t}{h^2} & 1+\frac{2b*\Delta t}{h^2}\\ &... &... &...\\ &&...&...&...&\\ &&&&&\\ &&&...&...&...&\\ &&&&...&...&...&\\ &&&&&-\frac{b*\Delta t}{h^2}&1+\frac{2b*\Delta t}{h^2} \end{bmatrix} A=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡1+h22b∗Δt−h2b∗Δt−h2b∗Δt1+h22b∗Δt−h2b∗Δt...−h2b∗Δt1+h22b∗Δt..............................−h2b∗Δt...1+h22b∗Δt⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤

令AU=C,其中：
U=(U10,U20,...UI−11)=(U^0_1,U^0_2,...U^1_{I-1})=(U10,U20,...UI−11)
C=(b∗Δth2∗U01+U10,U20,...UI−20+b∗Δth2∗UI−10)=(\frac{b*\Delta t}{h^2}*U^1_0+U^0_1,U^0_2,...U^0_{I-2}+\frac{b*\Delta t}{h^2}*U^0_{I-1})=(h2b∗Δt∗U01+U10,U20,...UI−20+h2b∗Δt∗UI−10)

可解出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);      endb = 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/2plot([0:h:1],sol(:,j),'r')axis([0 1 0 1]);xlabel('x');ylabel('u');hold on;   elseif j==Nplot([0:h:1],sol(:,j),'y')axis([0 1 0 1]);xlabel('x');ylabel('u');hold on;endlegend('t=0.001s','t=0.05s','t=0.1s');
end

在命令行输入：

得到结果矩阵：

得到二维图像：

画出三维图像：

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

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

隐方程求解一维抛物型方程（热传导方程）相关推荐

一维抛物型方程的差分解法
一维抛物型方程的差分解法问题内容算法求解显格式(explicit scheme)求解隐格式(implicit scheme)求解 Crank_Nicolsm格式求解实验结果三格式的迭代求解 ...
傅里叶谱方法-傅里叶谱方法求解一维 KdV 方程及其Matlab程序实现
3.3 傅里叶谱方法求解复杂偏微分方程 (组) 3.3.1 一维 KdV 方程背景介绍科特韦赫-德弗里斯方程(英语:Korteweg-De Vries equation),一般简称KdV方程,是1 ...
深度学习求解一维burgers方程和Galerkin求解泊松方程
一维burgers方程介绍 { u t + u u x − ( 0.01 / π ) u x x
打靶法c语言程序,打靶法求解一维Schrodinger方程程序示例
流程图 # include # include # defineh 10 floatksqu(float x,float e) {float k=0.1,m=0.01,ksq,v; v=k*(x-1) ...
【组合数学】递推方程 ( 常系数线性非齐次递推方程求解 | 递推方程标准型及通解 | 递推方程通解证明 )
文章目录一.递推方程标准型及通解二.递推方程通解证明一.递推方程标准型及通解 H(n)−a1H(n−1)−⋯−akH(n−k)=f(n)H(n) - a_1H(n-1) - \cdots - a ...
抛物型方程的差分解法matlab,抛物型方程的差分解法
? 0 时 2 为前差方程,当 ? ?1 时为后差方程.用控制体积法构造差分方程总是守恒型差分方程. (4) 积分方法采用积分方法构造差分方程基本思想是把微分...... 第7卷第4期2008年12 ...
有限元方法求解一维扩散方程(FEALPy)
有限元方法求解一维扩散方程文章目录有限元方法求解一维扩散方程有限元方法推导差分格式的介入数值算例之前完成了 FEALPy 有限元求解 Poisson 方程的数值算例, 通过湘潭大学王唯 ...
matlab riccati 方程,matlab解riccati方程
Riccati 方程求解 ? ? 前三个问题将介绍解析解与数值解,后一个属于非... 矩阵的化零空间或基础解系计算,支持符号运算求解连续 Lyapunov 方程.Sylvester 方程的数值解 ...
matlab将求解sin隐式解,Matlab隐式符号方程求解和赋值
近日处理了一个隐式方程的求解,由于方程含有较多的未知数,而且这些参数均是跟实验相关的一些参数,所以,必须得到需要求解的解与这些参数之间的一个表达式.之前是考虑用的Maple推导求解了该隐私方程,求解 ...

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

隐方程求解一维抛物型方程（热传导方程）相关推荐

最新文章

热门文章