抛物形扩散方程的有限差分法及数值实例
扩散方程的差分解法
（24）
则
， （25）
则误差方程为
（26）
误差放大因子为
（27）
要满足稳定性条件，则要求对所有的k值均有 。从（28）式中可以看出，当 （即 ）时， 恒成立。因此，全隐格式是无条件稳定的。
4.4收敛性
如果差分方程的解为 ，微分方程的解为 ，若当 ， 时，差分方程的解与微分方程的解之差
扩散方程的差分解法
在研究热传导过程、扩散过程、边界层现象时，我们常常遇到抛物型方程，这类方程中最典型、最简单的就是热传导方程。热传导方程中的自变量中包括时间t，它是描述一种随时间变化的物理过程，即所谓不定常现象。这类问题的基本定解问题应是初值问题，即在初始时刻（t=0）时给定定解条件，求解t>0时的解。
write(2,*) 'x=',x,'m'
do n=1,nt,200
write(2,*) (n-1)*dt,u(j,n)
enddo
if
enddo
!-----------!
end
5.3.2全隐格式
!----------------------------------------全隐格式求解扩散方程-----------------------------------------------!
由以上对一维扩散问题的分析，可知，求解一维扩散方程需给定初始条件及边界条件。
在本文计算中，取 ， 。
初始条件（ 时）
（29）
边界条件为
（30）
其初始时刻（ ）时的u分布如图1所示，x=0m处u随时间变化情况如图2所示，x=10m处u随时间变化情况如图3所示。
图1初始时刻u分布图
图2 x=0处u随时间变化图
抛物型方程的差分方法
抛物型方程的差分方法抛物型方程是描述物理现象中的薄膜振动、热传导、扩散等过程的方程，具有非常重要的应用价值。
差分方法是一种常用的数值计算方法，用于求解微分方程，对于抛物型方程的数值求解也是非常有效的方法之一、本文将介绍抛物型方程的差分方法，并具体讨论用差分方法求解抛物型方程的一些具体问题。
首先，我们来介绍一下抛物型方程的一般形式。
抛物型方程一般可以表示为：∂u/∂t=α(∂^2u/∂x^2+∂^2u/∂y^2)其中，u(x,y,t)是待求函数，t是时间，x和y是空间变量，α是常数。
这个方程描述的是物理过程中的扩散现象，如热传导过程、溶质的扩散过程等。
差分方法的基本思想是将求解区域离散化为一个个网格点，然后在每个网格点处用近似的方式来计算待求函数的值。
差分方法的求解步骤主要包括以下几个方面：1.选择适当的网格和步长。
在求解抛物型方程时，需要确定空间变量x和y所在的网格点以及步长，同时也需要确定时间变量t所在的网格点和步长。
通常，我们会选择均匀网格，步长选择合适的值。
2.建立差分格式。
差分格式是差分方法的核心部分，它包括对方程进行近似处理和离散化。
对于抛物型方程，常用的差分格式有显式差分格式和隐式差分格式等。
其中，显式差分格式的计算速度快，但是有一定的稳定性限制，而隐式差分格式的稳定性较好，但是计算量较大。
因此，在具体问题中需要根据实际情况选择适当的差分格式。
3.编写计算程序。
在建立差分格式后，需要编写计算代码来求解离散方程。
具体编写的过程包括定义初始条件、建立迭代计算过程、以及计算结果的输出等。
4.计算结果的验证与分析。
求解方程后，需要对计算结果进行验证和分析，主要包括对数值解和解析解的比较、对误差的估计和控制等。
在具体求解抛物型方程时，还会遇到一些问题，例如边界条件的处理、稳定性和收敛性的分析等。
下面将对其中一些问题进行详细讨论。
1.边界条件的处理。
边界条件对差分格式的求解结果有着重要的影响，常见的边界条件包括固定端（Dirichlet）边界条件和自由端（Neumann）边界条件等。
10_抛物型方程的有限差分方法
10_抛物型方程的有限差分方法抛物型方程是一类常见的偏微分方程，广泛应用于自然科学和工程学的领域中。
有限差分方法是一种常用的数值求解抛物型方程的方法之一、本文将介绍抛物型方程的有限差分方法（II）。
有限差分方法主要基于离散化的思想，将偏微分方程转化为差分方程，进而求解差分方程的数值解。
对于抛物型方程，其一般形式可以表示为：∂u/∂t=Δu+f(x,t)其中，u(x, t)是未知函数，表示空间位置x和时间t上的解，Δu表示Laplace算子作用于u的结果，f(x, t)是已知函数。
有限差分方法的基本思想是将空间和时间域进行离散化，将连续的空间和时间划分为有限个网格点，然后使用差分近似代替偏导数，得到差分方程。
假设空间域被划分为Nx个网格点，时间域被划分为Nt个网格点，对于每个网格点(i,j)，可以表示为(x_i,t_j)，其中i=0,1,...,Nx，j=0,1,...,Nt。
在有限差分方法中，我们使用中心差分近似来代替偏导数。
对于时间导数，可以使用向前差分或向后差分，这里我们使用向前差分，即：∂u/∂t≈(u_i,j+1-u_i,j)/Δt对于空间导数，可以使用中心差分，即：∂^2u/∂x^2≈(u_i-1,j-2u_i,j+u_i+1,j)/Δx^2将上述差分近似代入抛物型方程中，可以得到差分方程的离散形式：(u_i,j+1-u_i,j)/Δt=(u_i-1,j-2u_i,j+u_i+1,j)/Δx^2+f_i,j其中，f_i,j=f(x_i,t_j)。
重排上式，可以得到递推关系式：u_i,j+1=αu_i-1,j+(1-2α)u_i,j+αu_i+1,j+Δt*f_i,j其中，α=Δt/Δx^2通过设置初始条件和边界条件，可以利用以上递推关系式得到抛物型方程的数值解。
总结来说，抛物型方程的有限差分方法（II）是一种常用的数值求解抛物型方程的方法。
它基于离散化的思想，将偏微分方程转化为差分方程，然后利用中心差分近似代替偏导数，得到差分方程的离散形式。
(十二章)抛物型方程有限差分法
（3.3a） （3.3b）
； 其中
， LOD格式的计算步骤可以总结如下：
1） 令，。 2） 求解三对角线性方程组（3.3a）得到差分解。 3） 若，则增加1，转步骤4）。否则转4）。 4） 令。 5） 求解三对角线性方程组（3.3b）得到差分解。 6） 若，则增加1，转步骤5）。否则转7）。 7） 若，则增加1，转步骤2）。否则结束。
,
及边值条件
，
假定和在相应的区域光滑，并且于，两点满足相容条件,
则上述问题有唯一的充分光滑的解。
现在考虑边值问题（1.1），（1.3）的差分逼近
取 为空间步长，为时间步长，其中，是自然数，
,； ,
将矩形域分割成矩形网格。其中 表示网格节点；
表示网格内点（位于开矩形中的网格节点）的集合；
表示位于闭矩形中的网格节点的集合；
时我们简单地称差分格式稳定。
冯诺依曼稳定性分析估量了误差的放大或扩大。对一
种稳定的方法，必须选取步长使误差的放大因子不大于1.
前面讨论的向前差分格式（1.4）当网比时稳定，当时不稳定。这就
意味着给定空间步长以后，时间步长必须足够小，才能保证稳定。而向
后差分格式（1.6）和Grank-Nicholson格式（1.8）则对任何网比都是稳
+=+ ，==0。求出，在由取，可利用，解出，。如此下去，即可逐 层算出所有，。
如此每层必须解一个三对角线性方程组的格式称为隐格 式。并视为的近似值。
直观地说，采用显式格式进行求解既方便又省工作量。但 是，后面我们将看到，有些情况用隐式格式更为便利。
1.2.3 Grank-Nicholson法 将向前差分格式和向后差分格式做算术平均，得到的差分 格式称之为六点对称格式，也称为Grank-Nicholson格式：
抛物方程的有限差分法
图1
,我们需要求解这1/h +1()×T/τ+1()个点对应的函数值实上由已知的初边值条件蓝色标记附近的点可直接得到,所以只要确定微分方程的解在其它点上的取值即可,可记为u []
k j
=u (x j ,t k )。
建立差分格式
j =1, (1)
-1;k =0,1,…,T τ-1,用向前差分代替关于时间的
一阶偏导数,用二阶中心差分代替关于空间的二阶偏导数,则可定义最简显格式:
-u k j =u k j+1-2u k j +u k
j-1
h
2
变形有:
(上接第50页)极大值理论,检测初始行波、故障点反射波和对端母线反射波到达测量端的时间,测量故障点距离,从测试结果看,该方案有效弥补传统行波测距的不足之处,提高了故障测距的精确度。
【参考文献】
[1]陈靖.行波法故障测距的理论研究及其实现方案[D].武汉:武汉大学,2004.数值解的剖分图如图2:
图2
真解与数值解的误差剖分图如图3:
图3
3数值实验及结果分析
我们对所求解的初边值问题(1)进行算法精度的数值实验,当
u 0
(x )sin πx 时,边界值仍然为u (0,t )=u (1,t )=0,其精确解为:u (x ,t )
从表中我们可以看出。
. All Rights Reserved.。
抛物型方程的差分方法
2) 逆风差分格式
1 un un j j
a
n un j u j 1
§5 抛物型方程的差分方法
本章內容: 1． 常系数抛物型方程的初值问题 2．初边值问题 3．对流扩散方程 4．Richardson外推法 ^^数值实验4（网格比的重要性） ^^数值实验5 (Richardson 外推的精度) 5．变系数方程 6．二维抛物型方程问题的计算
1． 常系数抛物型方程初值问题
1 2
返回本节
c. 三层显式关系 Richardson格式 n n 1 1 un un un j j j 1 2u j u j 1 a 0 h2 2 Du Fort－Frankel格式
1 1 un un j j
2
a
n 1 1 n un un j 1 (u j j ) u j 1
目标点：Jh 1, x j jh, N T , ( x j , t n )
1 n n un un un j j j 1 2u j u j 1 内部点的离散： a 0 h2
边值点的离散： u0 ( n ), u J (n )
n n
初值点的离散： u j f (x j ) f j
1 4(1 )a sin 2 1 4a sin 2 kh 2
kh 2
3
G ( , k ) 1 1 G ( , k ) 1 4a (1 2 )sin 2 2a (1 2 ) 1 kh 2 2 0 1 2
有限差分求解扩散的数值解
Jeff
C2
无孔结构
J0 =-D ∂C ∂z
多孔介质
ε ∂C Jeff =- D τ ∂z
其中ε为孔隙率
曲折因子
τ= J0 ε Jeff
3.扩散模拟表征多孔介质物质输运性能（相转换流延制备） SRCT构建孔结构数组
灰度图
二值图
黑色：孔 白色：实体 孔隙率：0.35
三维结构重构图A
3.扩散模拟表征多孔介质物质输运性能
1.有限差分的基本原理（ 一维扩散例子） 计算结果
稳态下
∂C ∂2 C =D 2 = 0 ∂t ∂x
2.边界条件处理（ 二维柱坐标扩散例子）
扩散进口
扩散出口
绝缘 环
圆柱上底面和侧面为高浓度面（红色） 圆柱下底面同心小圆为低浓度面（白色） 圆柱下底面圆环为绝缘面 （蓝色）
2.边界条件处理（ 二维柱坐标扩散例子） 0 Z
r
扩散微分方程： ∂C ∂2 C 1 ∂C ∂2 C =D( 2 + r + 2) ∂t ∂r ∂r ∂Z
圆心对称
柱坐标下取体积微元由质量衡算推导扩散偏微分方程
r+ △r
△Z
r
体积微元质量变化量
∆m=-2π r+∆r ∆Z∆t*J -π r+∆r 2 -r2 ∆t*J
移项化简得
r+∆r+2πr∆Z∆t*J r Z
2 2 + π r+ ∆r -r ∆t*J Z+∆Z
∆m 1 -J = * 2πr∆Z∆r∆z r
r+∆r *
即
r+∆r -J r *r J Z+∆Z -J Z ∆r ∆Z ∂C ，将 J=-D 代入可得 ∂x
2.2 抛物型方程的差分解法
u ( j 1, n) 2u ( j , n) u ( j 1, n) u 2h 4 ( j , n) u ( j , n) 2 2 4 h 2 t 4! x
n
（8）
0
Lu j
n
Lh, u j R j n
式中：
2 4 2 2 h 2 Rn u ( j , n ) u ( j , n ) O ( h ) j 4 2 4! x 2 t
(backward space difference) (backward time difference)
u n j
（3）一阶中心差分（central difference）
hu
n j
un 1 un
j 2
j
1 2
h
u
n j
uj
n
1 2
uj
n
1 2
1 n 1 un u j j
n
（22）
n+1 n
j-1
j
j+1
注意：
① 泰勒展开点在格边上，不是在结点上，但在格式中未出现格边量。 ② ③
O( 2 h2 ) ——全二阶精度。 1 在 ( j, n ) 点展开时，用到了周围6个结点上的量，该格式又称为六点格式。 2 Rj
2u idea：是将微分方程中的 2 项以 u ( x, t ) x
u j n1 u j n 1 2
u j 1n 2u j n u j -1n h2
0 (23)
(24)
u j n1 2r(u j 1n - 2u j n u j 1n ) u j n1
扩散模型的原理范文
扩散模型的原理范文扩散模型是一种模拟物质传播和扩散过程的数学模型，广泛应用于社会科学、物理学、生物学、工程学等领域。
其基本原理是描述和预测物质在空间和时间上的分布变化规律。
以下将详细介绍扩散模型的原理。
1.扩散过程扩散是指物质由高浓度区域向低浓度区域的自发传播过程。
这种传播可以是由分子、粒子、信息等传递激发的，也可以是由温度、浓度梯度等驱动的。
扩散过程遵循费克定律，即物质扩散的速率与浓度梯度成正比。
2.扩散方程∂C/∂t=D∂²C/∂x²其中，C是物质浓度，t是时间，x是空间位置，D是扩散系数。
该方程描述了物质浓度随时间变化的速率和空间扩散的强度。
对于多维情况，扩散方程的形式会有所不同。
3.初始条件与边界条件扩散模型还需要指定初始条件和边界条件。
初始条件是指在初始时刻，物质浓度在空间上的分布情况。
边界条件则是指在边界位置上，物质浓度的变化规律。
这些条件通常通过实验测量或数值模拟获得。
4.数值求解由于扩散方程是偏微分方程，对于大多数情况，无法求得解析解。
因此，需要使用数值方法来近似求解扩散方程。
常见的数值方法包括有限差分法、有限元法等。
这些方法将空间和时间分割成小的网格，通过在网格上进行计算来模拟物质扩散过程。
5.扩散系数和扩散机制扩散系数是一个重要的参数，它反映了物质在扩散过程中的扩散速率。
扩散系数的大小取决于物质本身的性质，以及扩散过程中的环境条件。
不同的扩散机制也会导致不同的扩散系数。
例如，在气体中，扩散主要由分子碰撞引起；在液体中，扩散主要由颗粒之间的扩散引起。
6.扩散模型的应用扩散模型在实际中有广泛的应用。
在社会科学中，扩散模型可以用于描述和预测信息、流行病等的传播过程。
在物理学中，扩散模型可以用于研究热传导、原子扩散等现象。
在生物学中，扩散模型可以用于模拟细胞内物质的运输过程。
在工程学中，扩散模型可以用于设计和优化化学反应器、催化剂等。
总结：扩散模型是一种用于模拟物质传播和扩散过程的数学模型。
第四章 抛物型方程的有限差分方法
2 h 称为Du Fort -Frankel格式，仍为三层显式格式．
2
a
n 1 n 1 n un ( u u ) u j 1 j j j 1
0
截断误差： T x j , tn a u x j , tn u x j , tn 2 u x j h, tn u x j , tn u x j , tn u x j h, tn h2
1 2a G , k 0
0 4a cos kh 1 2a 1 1 0 4a cos kh 1 2a 1 2a 1 2a 0 1
2
1
4a cos kh 2a 1 G , k 的特征方程： 0 1 2a 1 2a
修正 Richardson：无条件不稳定显格式
Du Fort Frankel：无条件稳定的三层显格式． 但后者的相容性是有条件的．事实上， 显格式中，无条件相容和无条件稳定是无法同时成立的．
4 三层隐式格式
先考虑
n 1 n u u 3 j j n n 1 u u 1 j j 1 n1 n1 un 2 u u j 1 j j 1
引理1.1实系数二次方程 2 b c 0的根： c 1． 模 1 b 1 c， " "设1 , 2是方程两根，且 i 1 i 1, 2 证： c b 则12 c1 2 b a a 12 c c 1 2 1 1 2 b 1 c b 1 12 1 2 1 12 1 2 1 1 1 2 0, 若 1 2 0 1 12 1 2 1 1 1 2 0, 若1 2 0 b 1 c
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
偏微分方程数值解所在学院：数学与统计学院课题名称：抛物形扩散方程的有限差分法及数值实例学生姓名：**抛物形扩散方程的有限差分法及数值实例1.1抛物型扩散方程抛物型偏微分方程是一类重要的偏微分方程。
考虑一维热传导方程：22(),0u ua f x t T t x∂∂=+<≤∂∂ (1.1.1) 其中a 是常数，()f x 是给定的连续函数。
按照初边值条件的不同给法，可将(1.1.1)的定解分为两类：第一，初值问题（Cauchy 问题）：求足够光滑的函数()t x u ,，满足方程(1.1.1)和初始条件：()()x x u ϕ=0,, ∞<<∞-x (1.1.2)第二，初边值问题（也称混合问题）：求足够光滑的函数()t x u ,，满足方程(1.1.1)和初始条件：()()x x u ϕ=0,, 0x l << (1.1.3) 及边值条件()()0,,0==t l u t u ， T t ≤≤0 (1.1.4)假定()x f 和()x ϕ在相应的区域光滑，并且于()0,0，()0,l 两点满足相容条件,则上述问题有唯一的充分光滑的解。
1.2抛物线扩散方程的求解下面考虑如下热传导方程22()(0.)(,)0(,0)()u ua f x t x u t u L t u x x ϕ⎧∂∂=+⎪∂∂⎪⎪==⎨⎪=⎪⎪⎩(1.2.1) 其中，0x l <<，T t ≤≤0,a (常数)是扩散系数。
取N l h =为空间步长，MT=τ为时间步长，其中N ，M 是自然数，用两族平行直线jh x x j ==,()N j ,,1,0 =和k t t k τ==, ()M k ,,1,0 =将矩形域G {}T t l x ≤≤≤≤=0;0分割成矩形网格。
其中 (),j k x t 表示网格节点；h G 表示网格内点（位于开矩形G 中的网格节点）的集合；h G 表示位于闭矩形G 中的网格节点的集合；h Γ表示h G -h G 网格边界点的集合。
k j u 表示定义在网点(),j k x t 处的待求近似解，N j ≤≤0，M k ≤≤0。
现在对方程进行差分近似： (一) 向前差分格式=-+τk jk j u u 11122(())k k kj j j j j j u u u af f f x h+--++= (1.2.2)()j j j x u ϕϕ==0， k u 0=k N u =0 (1.2.3)计算后得：111(12)k k k k j j j j j u ru r u ru f τ++-=+-++ (1.2.4)其中，2a r h τ=，1,,1,0-=N j ，1,,1,0-=M k 。
显然，这是一个四点显示格式，每一层各个节点上的值是通过一个方程组求解到的。
方程组如下：10001210110002321210003432310001121(12)(12)(12)(12)N N N N N u ru r u ru f u ru r u ru f u ru r u ru f u ru r u ru f ττττ----⎧=+-++⎪=+-++⎪⎪=+-++⎨⎪⎪⎪=+-++⎩ (1.2.5) 若记()TkN k k k u u u 121,,,-= u ，()()()()T N x x x 121,,,-=ϕϕϕϕ ，()()()()T N x f x f x f 121,,,-=τττ f 则显格式(1.2.4)可写成向量形式10,0,1,,1k k k M φ+⎧=+=-⎨=⎩u Au f u (1.2.6)其中⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛----=r r r r r r r rr r21002100210021A而对于向前差分格式，当网比12r ≤时稳定，当12r >时不稳定。
这就意味着给定空间步长h 以后，时间步长τ必须足够小，才能保证稳定。
1.3抛物型热传导方程数值算例对于(1.2.1)所描述的扩散方程，取1a =已知方程的精确解为2sin t u e x ππ-= ：22(,0)sin()(0,)(1,)001,00.5u ut x u x x u t u t x t π⎧∂∂=⎪∂∂⎪⎪=⎨⎪==⎪≤≤≤≤⎪⎩ (1.3.1) 设空间步长1/h M =，时间步长为0.5/N τ=，网格比2/r g h =。
向前格式：11122,1,...,1,1,...,k k k k kj jj j j u u u u u j M k N hτ++---+==-=边值条件：110010111220,,k kk k kk k u u u u u j u u h τ++----+===， 11111122,,k k k k k k k M MM M M M M u u u u u j M u u hτ+++-+---+===. 初值条件：(,0)sin ,0,1,...,j j u x x j M π==对时间和空间进行分割，令M=40，N=1600，通过Matlab 计算得到该方程的解析解，数值解以及相对误差如下：图(1)解析解的图像图(2)数值解的图像图(3)M=40,N=1000的相对误差的图像我们取部分精确解和数值解进行比较，结果如表(1)x t数值解精确解相对误差0.1 0.1 0.11510.115241.170010-⨯0.2 0.2 0.08150.08164⨯1.658010-0.3 0.3 0.04180.04194⨯1.275210-0.4 0.4 0.01830.01845⨯7.445610-0.5 0.45 0.01170.01185⨯5.375510-0.6 0.35 0.03000.030141.067410-⨯0.7 0.25 0.06840.06864⨯1.741010-0.8 0.15 0.50840.50855⨯7.589710-0.9 0.05 0.36540.36545⨯1.740710-表(1)数值解与精确解的比较由表(1)我们可以看出，精确解和数值解的绝对误差在410-以内，因此可以得出，在分割M=40，N=1600下，该有限差分方法对方程(1.3.1)是收敛和稳定的。
下面，我们比较在不同的分割下对有限差分算法精度的影响。
在扩散系数1a=不变的情况下，讲时间和空间进行更加细密的分割，取==，其中，M表示空间上的分割，N表示时间上的分割。
观察数M N50,10000值解与精确解在节点,x t处的绝对误差值，如下图所示：j k图(4)M=50,N=10000的相对误差的图像由图(3)和图(4)，两者在节点处的误差收敛分别是在410-和510-以内，因此，可以得出的结论是：在收敛范围内，随着时间和空间的分割越细，节点数越多，精确解和解析解之间的绝对误差也越小，有限差分法的算法精度也越高。
最后，我们比较网比1/2r=以及1r=时扩散方程的收敛情况。
当网比1r=时，此时我们取M=10，N=50，这时，方程的数值解与解析解还有相对误差图如下：图(5)M=10,N=50的解析解的图像图(6)M=10,N=50的数值解的图像图(7)M=10,N=50的绝对误差的图像此时，我们观察绝对误差发现，扩散方程(1.3.1)时不收敛不稳定的。
而前面我们已经知道，到网格比为12r=时，方程是收敛稳定的。
所以，我们可以验证，当网比12r≤时稳定，当12r>时不稳定。
[参考文献][1]李荣华，刘播.微分方程数值解法[M].北京.高等教育出版社.2009.1.[2]王曰朋.偏微分方程数值解[OL]./link?url=-UAAXRI8_V547zVS6UwenT8rGKsbwNf9CuDmh2qmsy5K8eQ 32PzhYazZ_sWiHfz_Pj3LA7ufHH4O3t8NIlUCnXNUiyYhssqmNBeArsrLwQG[3]未知.偏微分方程的Matlab解法[OL]./link?url=WBlR0q9n02YsNg6UOLkCQ4Z5QOCefDKNdEglrNR2TJ8Vlxa J9IkCrLlaEDlJ_OHCLehf19UJU_B7PRe1skQTYaaLcAa1UWcwlv7GYsWNtG7[4]周品，何正风.MATLAB数值分析.[M].北京.机械工业出版社.2009.1.附录：L=1;M=40;N=1600;alfa=1;lambda=0.5;%网格比%**********************************************%h=L/M;%空间步长x=0:h:L;x=x';tao=lambda*h^2/alfa;%时间步长tm=N*tao;%热传导的总时间tm%tm=0.1;t=0:tao:tm;t=t';%计算初值和边值U=zeros(M+1,N+1);U(:,1)=sin(pi*x);U(1,:)=0;U(M+1,:)=0;%*************用差分法求出温度U，与杆长L，时间T的关系****************% for k=1:Nj=2;while j<=MU(j,k+1)=lambda*U(j+1,k)+(1-2*lambda)*U(j,k)+lambda*U(j-1,k);j=j+1;endendlength(U);%***************设置立体网格****************%for i=1:N+1X(:,i)=x;endfor j=1:M+1Y(j,:)=t;endmesh(X,Y,U);legend('数值解');xlabel('X');ylabel('T');zlabel('U');z=zeros(M+1,N+1);for j=1:M+1for k=1:N+1z(j,k)=exp(-pi*pi*t(k))*sin(pi*x(j));endend%mesh(x,t,z')legend('解析解');xlabel('X');ylabel('T');zlabel('Z');for j=1:M+1for k=1:N+1y(j,k)=abs(z(j,k)-U(j,k));endendmesh(x,t,y');legend('绝对误差');xlabel('X');ylabel('Y');zlabel('error');。