一维对流扩散方程的数值解法
- 格式:pdf
- 大小:374.54 KB
- 文档页数:6
一维非稳态对流-扩散方程的隐式中心差分格式一维非稳态对流-扩散方程描述了在一维空间中同时存在对流和扩散过程的物理现象。
这个方程在很多工程和科学领域都有广泛的应用,如热传导、质量传输和流体力学等。
对方程进行数值求解可以得到物理现象的定量解,进而对系统行为进行预测和优化。
对于一维非稳态对流-扩散方程的数值解法中,中心差分格式是一种常用的方法。
中心差分格式是基于中心差分近似的方法,该方法精确地处理了对流和扩散效应,并适用于广泛的问题。
其中,隐式格式是一种特殊的中心差分格式,它在处理高度非稳态情况下的数值解具有优势。
在一维非稳态对流-扩散方程的隐式中心差分格式中,我们假设空间网格的节点为x_i,时间步长为Δt。
方程的解u(x,t)在网格节点x_i处的近似值为u_i^n,其中n表示时间步数。
根据对流和扩散项的中心差分近似,方程可以离散为如下的格式:(u_i^{n+1}-u_i^n)/Δt=α(u_{i-1}^{n+1}-2u_i^{n+1}+u_{i+1}^{n+1})/(Δx^2)-β(u_{i+1}^{n+1}-u_{i-1}^{n+1})/(2Δx)其中,α表示扩散系数,β表示对流系数,Δx表示空间步长。
在隐式中心差分格式中,时间步n+1的解u_i^{n+1}是未知的,我们将其视为待求解的值。
通过将方程的右侧扩散和对流项全部取为n+1步的值,从而得到一个关于u_i^{n+1}的线性方程。
因此,我们可以得到如下的表达式:u_i^{n+1}=(αΔt/(Δx^2)-βΔt/(2Δx))u_{i-1}^{n+1}+(1-2αΔt/(Δx^2))u_i^{n+1}+(αΔt/(Δx^2)+βΔt/(2Δx))u_{i+1}^{n+ 1}这个方程可以用矩阵的形式表示为:AU^{n+1}=BU^n其中,U^{n+1}是一个列向量,包含了所有网格节点i处的解u_i^{n+1};U^n是一个列向量,包含了所有网格节点i处的解u_i^n;A和B是相关系数矩阵,具体的表达式为:A=[1-2αΔt/(Δx^2),αΔt/(Δx^2)+βΔt/(2Δx),0, 0[αΔt/(Δx^2)-βΔt/(2Δx),1-2αΔt/(Δx^2),αΔt/(Δx^2)+βΔt/(2Δx), 0[0,αΔt/(Δx^2)-βΔt/(2Δx),1-2αΔt/(Δx^2), 0...[0,...,0,αΔt/(Δx^2)-βΔt/(2Δx),1-2αΔt/(Δx^2)]B = identity matrix通过求解这个线性方程组,就可以得到隐式中心差分格式下的数值解。
一维对流扩散方程是指一维均匀的边界层上的传质过程的数学模型,常用于描述对流扩散过程中的温度、湿度、速度等场的分布情况。
一维对流扩散方程的数学形式为:∂φ/∂t+U∂φ/∂x=D∂^2φ/∂x^2其中φ表示传质物质的浓度,t表示时间,x表示空间坐标,U表示对流速度,D表示扩散系数。
二维对流扩散方程是指二维均匀的边界层上的传质过程的数学模型,常用于描述对流扩散过程中的温度、湿度、速度等场的分布情况。
二维对流扩散方程的数学形式为:∂φ/∂t+U∂φ/∂x+V∂φ/∂y=D∂^2φ/∂x^2+D∂^2φ/∂y^2其中φ表示传质物质的浓度,t表示时间,x和y分别表示两个空间坐标,U和V分别表示两个方向上的对流速度,D表示扩散系数。
单调差分格式是一种常用的数值求解方法,它通过进行差分运算来求解微分方程的数值解。
在求解一维和二维对流扩散方程时,可以使用单调差分格式来解决。
具体来说,可以将空间坐标和时间分别离散化,将对流扩散方程转化为一个线性方程组,然后使用单调差分格式来解决。
单调差分格式的具体形式取决于方程的类型和离散化的方式,但一般来说,它都是将微分方程的差分形式写成一个线性方程组的形式。
例如,在求解一维对流扩散方程时,可以使用下面的单调差分格式:φ_i^{n+1}=φ_i^n+Δt(D(φ_{i+1}^n-2φ_i^n+φ_{i-1}^n)/Δx^2+U(φ_ {i+1}^n-φ_{i-1}^n)/2Δx)其中φ_i^n表示第i个网格点在时间步n的浓度值,Δx和Δt分别表示网格的空间步长和时间步长。
同样的,在求解二维对流扩散方程时,可以使用下面的单调差分格式:φ_i^n=φ_i^n+Δt(D(φ_{i+1,j}^n+φ_{i-1,j}^n+φ_{i,j+1}^n+φ_{i,j-1}^ n-4φ_i^n)/Δx^2+U(φ_{i+1,j}^n-φ_{i-1,j}^n)/2Δx+V(φ_{i,j+1}^n-φ_ {i,j-1}^n)/2Δy)其中φ_i^n表示第(i,j)个网格点在时间步n的浓度值,Δx和Δy分别表示网格在x和y方向上的空间步长,Δt表示时间步长。
对流方程及其解法对流方程是描述流体运动的最基本方程之一,涉及热、动量、物质等的传递现象,对于各种物理问题的研究都具有重要意义。
本文将从对流方程的基本形式和意义出发,探讨其常见解法及相关应用。
一、对流方程的基本形式与意义对流方程是描述流体中质量、热量和动量传递的方程,其基本形式可以写作:$$ \frac{\partial\phi}{\partial t} + (\mathbf{v}\cdot\nabla)\phi =\nabla\cdot(\Gamma\nabla\phi) $$其中,$\phi$为描述流体量的变量,如温度、密度、浓度等;$\mathbf{v}$为流体的流速,$\Gamma$为扩散系数。
对该方程的解析求解较为困难,故通常采用数值方法进行求解。
下面介绍几种常见的数值解法。
二、有限差分法有限差分法是在连续方程的基础上,利用有限差分代替导数,将微分方程变为代数方程组,从而利用计算机求解的方法。
其基本思想是将求解区域划分为有限个网格,对每个网格内的量用差分代替导数,从而得到有限差分方程。
以简单的二维对流扩散为例,其对流方程为:$$ \frac{\partial\phi}{\partial t} + u\frac{\partial\phi}{\partial x} + v\frac{\partial\phi}{\partial y} = \Gamma\frac{\partial^2\phi}{\partial x^2} + \Gamma\frac{\partial^2\phi}{\partial y^2} $$其中,$u$和$v$分别代表$x$和$y$方向的流速。
对该方程进行离散,假设$\phi_{i,j}$为$x=i\Delta x$,$y=j\Delta y$处的$\phi$值,则可以得到:$$ \frac{\phi^{k+1}_{i,j} - \phi^k_{i,j}}{\Delta t} +u\frac{\phi^k_{i+1,j} - \phi^k_{i-1,j}}{2\Delta x} +v\frac{\phi^k_{i,j+1} - \phi^k_{i,j-1}}{2\Delta y} $$$$ = \frac{\Gamma\Delta t}{(\Delta x)^2}(\phi^k_{i+1,j} -2\phi^k_{i,j} + \phi^k_{i-1,j}) + \frac{\Gamma\Delta t}{(\Deltay)^2}(\phi^k_{i,j+1} - 2\phi^k_{i,j} + \phi^k_{i,j-1}) $$其中,$k$为时刻,$\Delta x$和$\Delta y$分别为$x$和$y$方向的网格间距。
Sinc-Chebyshev配置方法求解一维对流扩散方程毛志【摘要】利用复合移位Sinc函数和移位Chebyshev多项式,构造了求解变系数的一维对流扩散方程初边值问题的Sinc-Chebyshev配置方法.【期刊名称】《铜仁学院学报》【年(卷),期】2013(015)005【总页数】4页(P146-149)【关键词】Sinc函数;移位Chebyshev多项式;Sinc-Chebyshev配置方法;变系数的对流扩散方程【作者】毛志【作者单位】铜仁学院数学与计算机科学系,贵州铜仁554300【正文语种】中文【中图分类】O241.82一、引言对流扩散方程(convection diffusion equation)是一类基本的运动方程,是描述粘性流体的非线性方程的线性化模型方程。
它可以用来描述空气动力学、水力学、环境保护和生物、化学工程等众多科技和工程领域中的对流扩散问题[1],所以关于对流扩散方程数值方法的研究具有十分重要的理论价值和现实意义。
对流扩散问题的有效数值解法一直是计算数学中重要的研究内容。
由于对流扩散方程同时含有对流项和扩散项,在数值求解时会引起数值振荡和数值弥散[2],使得方程的求解比较复杂。
目前,求解对流扩散方程的数值方法有多种,如有限差分法(FDM)[3]、有限元法(FEM)[4][5]、有限体积法(FVM)[6][7]、边界元法(BEM)等,其中有限差分方法是一种重要的数值计算方法。
目前对于常系数的对流扩散方程已有较多、较好的研究,而对于变系数的对流扩散方程的研究却较少,需要对方程做一定的假设。
但在实际应用中,比如耦合流场的对流扩散方程,则需要研究变系数的情形[8]。
基于此,本文考虑如下变系数的一维对流扩散方程的初边值问题,其中a(x,t)、b(x,t)、f(x,t)和g(x)已知,a(x,t )≠0、 b (x,t )≠ 0 且连续。
本文设计了一种新的求解上述问题(1)~(3)的数值方法---Sinc-Chebyshev配置方法。
对流扩散方程有限差分方法流扩散方程是描述流体内部物质的扩散过程的方程,它可以用于描述溶质的扩散、热量的传导以及动量的传递。
在许多工程和科学领域中,比如地球科学、生物医学和工程学等,流扩散方程都有着广泛的应用。
在数值计算中,有限差分方法是一种常用的数值解法,可以非常有效地解决流扩散方程。
下面将详细介绍对流扩散方程有限差分方法的原理和步骤。
首先,考虑一维流扩散方程的一般形式:∂C/∂t=D∂²C/∂x²-V∂C/∂x其中,C是扩散物质的浓度,t是时间,x是空间位置,D是扩散系数,V是对流速度。
为了使用有限差分方法求解上述方程,我们需要将时间和空间分布离散化,得到方程在网格点上的近似表示。
首先,将时间轴分为n个等间隔的时间步长Δt,空间轴分为m个等间隔的网格点,网格点之间的间距为Δx。
然后,我们使用数值方法来逼近方程中的各个导数项,采用中心差分公式:∂C/∂t≈(C_i^(n+1)-C_i^n)/Δt∂²C/∂x²≈(C_i+1^n-2C_i^n+C_i-1^n)/Δx²∂C/∂x≈(C_i+1^n-C_i-1^n)/(2Δx)将上述近似代入流扩散方程,可以得到:(C_i^(n+1)-C_i^n)/Δt=D(C_i+1^n-2C_i^n+C_i-1^n)/Δx²-V(C_i+1^n-C_i-1^n)/(2Δx)整理上式,可以得到对流扩散方程的有限差分方程:C_i^(n+1)=C_i^n+(DΔt/Δx²)(C_i+1^n-2C_i^n+C_i-1^n)-(VΔt/2Δx)(C_i+1^n-C_i-1^n)上述方程给出了方程在时刻n+1时刻网格点i的值,即C_i^(n+1),它的值通过已知时刻n时刻各个网格点的值C_i^n来计算。
最后,我们可以使用迭代的方法,从初始条件C_i^0开始,依次计算下一个时刻的网格点C_i^(n+1),直到达到所需的计算精度或者计算到需要的时间步长。
一维扩散方程数值求解一维扩散方程是描述物质扩散过程的数学模型,广泛应用于物理、化学、生物和工程等领域。
本文将介绍一维扩散方程的数值求解方法,并探讨其在实际问题中的应用。
一维扩散方程的数值求解是通过离散化连续物理问题,将其转化为有限个代数方程的求解过程。
首先,我们需要将一维空间进行离散化,将其划分为一系列离散节点。
然后,通过数值方法近似计算节点上的物理量,如浓度、温度等。
最常用的数值方法包括有限差分法和有限元法。
有限差分法是一种简单且常用的数值求解方法。
它通过将偏导数用差商近似表示,将一维扩散方程转化为离散的代数方程组。
具体而言,我们可以使用向前差分、向后差分或中心差分等方式来近似计算偏导数。
然后,通过代数方程组的求解,得到离散节点上的物理量。
有限元法是一种更为灵活和精确的数值求解方法。
它将一维空间划分为一系列小单元,通过定义适当的插值函数,将节点上的物理量表示为有限个自由度的线性组合。
然后,通过求解线性方程组,得到每个单元上的物理量。
最后,通过汇总所有单元的解,得到整个一维空间上的物理量分布。
一维扩散方程的数值求解在许多领域都有广泛的应用。
在物理学中,它可以用于描述热传导、质量传递等过程。
在化学工程中,它可以用于模拟反应器内物质的传输与转化。
在生物学中,它可以用于研究细胞内物质的扩散行为。
在工程学中,它可以用于设计材料的扩散性能和优化结构。
除了基本的一维扩散方程,还可以考虑一些扩展问题。
例如,考虑非线性扩散系数、吸附效应、反应等因素。
这些扩展模型可以更准确地描述实际问题,但也增加了数值求解的难度。
一维扩散方程的数值求解是解决物质扩散问题的重要手段。
通过合理选择数值方法和适当的离散化方式,可以得到准确的物理量分布。
这为我们研究和解决实际问题提供了有力的工具。
同时,我们也需要注意数值误差和收敛性等问题,以确保数值结果的可靠性和有效性。
因此,深入理解一维扩散方程的数值求解方法,对于科学研究和工程应用都是非常重要的。
一维扩散方程差分格式的数值计算一维扩散方程是描述物质在一维空间中扩散过程的方程。
数值计算是一种近似求解微分方程的方法,可以通过离散化空间和时间来求解一维扩散方程。
本文将介绍一维扩散方程差分格式的数值计算方法,并给出一个具体的数值计算实例。
∂u/∂t=D∂²u/∂x²其中,u是扩散物质的浓度,t是时间,x是空间坐标,D是扩散系数。
差分格式的基本思想是将连续的时间和空间变量离散化为一系列有限的点,然后用离散化后的点代替原方程中的连续变量,从而得到一个差分方程。
一维扩散方程的差分格式数值计算方法有很多种,下面介绍两种基本的差分格式:显式差分格式和隐式差分格式。
1.显式差分格式:显式差分格式的基本思路是使用当前时间步的解来计算下一个时间步的解。
通过对一维扩散方程进行差分得到:(u_i)_(n+1)=(u_i)_n+D*(∆t/∆x²)*((u_(i-1))_n-2(u_i)_n+(u_(i+1))_n)其中,(u_i)_(n+1)表示时间步n+1时刻、位置i处的扩散物质浓度。
该公式使用当前时间步n的解来逐点计算下一个时间步n+1的解。
2.隐式差分格式:隐式差分格式的基本思路是使用下一个时间步的解来计算当前时间步的解。
通过对一维扩散方程进行差分得到:((u_i)_(n+1)-(u_i)_n)/∆t=D*(∆x²)*((u_(i-1))_(n+1)-2(u_i)_(n+1)+(u_(i+1))_(n+1))这是一个关于时间步n+1的隐式方程,需要使用迭代方法求解。
数值计算的实例:假设在一根长为L的杆上有一种扩散物质,杆的两端固定浓度为0,即u(0, t) = u(L, t) = 0;初始时刻杆上的浓度分布为一个正弦函数,即u(x, 0) = sin(πx/L);扩散系数为D。
我们需要计算杆上扩散物质的浓度随时间的变化情况。
首先,选择合适的网格间距∆x和时间步长∆t。
然后将杆上的空间坐标和时间离散化为一系列点,得到网格。
一维扩散方程差分格式的数值计算∂u/∂t=D∂²u/∂x²其中,u(x,t)是在位置x和时间t的扩散现象的浓度,D是扩散系数。
为了对一维扩散方程进行数值计算,可以使用差分格式。
最常用的差分格式是向前差分和中心差分。
1.向前差分格式:使用向前差分格式将时间t和位置x分别离散化,差分步长分别为Δt和Δx。
将扩散方程中的偏导数用有限差分近似替代,可以得到近似方程:(u_i(t+Δt)-u_i(t))/Δt=D(u_i-1(t)-2u_i(t)+u_i+1(t))/Δx²其中,u_i(t)表示在位置x_i和时间t的解,u_i(t+Δt)和u_i(t)是上一时刻和当前时刻的浓度,u_i-1(t)和u_i+1(t)分别是x_i左右两侧位置的解。
这样,一维扩散方程就被转化为一个差分方程。
根据初始条件u(x,0)和边界条件u(0,t)和u(L,t),L表示空间区域的长度,可以得到差分方程的初始条件。
使用向前差分格式可以得到一个显式迭代公式:u_i(t+Δt)=u_i(t)+DΔt(u_i-1(t)-2u_i(t)+u_i+1(t))/Δx²这个公式可以用来逐步推进时间t的步骤,从而获得扩散过程中的浓度分布。
2.中心差分格式:使用中心差分格式将时间t和位置x分别离散化,差分步长分别为Δt和Δx。
将扩散方程中的偏导数用有限差分近似替代,可以得到近似方程:(u_i(t+Δt)-u_i(t))/Δt=D(u_i-1(t)-2u_i(t)+u_i+1(t))/Δx²与向前差分格式不同的是,在右侧位置x_i+1处使用u_i+1(t)近似。
这个差分方程可以进一步简化为一个稳定的隐式迭代公式:u_i(t+Δt)=u_i(t)+DΔt(u_i-1(t+Δt)-2u_i(t+Δt)+u_i+1(t+Δt))/Δx²这个公式可以通过求解线性方程组来计算下一个时间步长的解。
以上是一维扩散方程差分格式的数值计算的基本原理和方法。
一维对流扩散方程是描述物质传输和扩散现象的重要数学模型,对于工程、地质、生物等领域具有重要的理论和应用价值。
在科学研究和工程实践中,人们经常需要利用计算机软件对一维对流扩散方程进行数值求解,以获得物质传输和扩散的详细信息。
MATLAB作为一种强大的科学计算软件,提供了丰富的数学工具和编程接口,可以方便地对一维对流扩散方程进行数值求解。
本文将介绍利用MATLAB对一维对流扩散方程进行数值求解的基本方法和步骤。
一、一维对流扩散方程的数学模型一维对流扩散方程是描述物质在一维空间中传输和扩散的数学模型,通常可以写成如下的形式:∂c/∂t + u∂c/∂x = D∂^2c/∂x^2其中,c是物质浓度,t是时间,x是空间坐标,u是对流速度,D是扩散系数。
该方程的求解可以得到物质浓度随时间和空间的变化规律,对于理解物质传输和扩散过程具有重要意义。
二、MATLAB求解一维对流扩散方程的基本步骤在MATLAB中,可以利用偏微分方程求解工具箱(Partial Differential Equation Toolbox)来对一维对流扩散方程进行数值求解。
求解的基本步骤如下:1. 网格的生成首先需要在空间上生成一个网格,将一维空间离散化为有限个网格点。
可以利用MATLAB中的linspace函数或者自定义函数来实现网格的生成。
2. 边界条件和初始条件的设定根据具体问题的边界条件和初始条件,需要在MATLAB中对边界条件和初始条件进行设定。
3. 偏微分方程的建立利用MATLAB中的偏微分方程建立工具箱,可以方便地将一维对流扩散方程建立为MATLAB中的偏微分方程对象。
4. 方程的数值求解利用MATLAB中的求解器对建立的偏微分方程进行数值求解,可以获得一维对流扩散方程的数值解。
5. 结果的可视化可以利用MATLAB中丰富的绘图函数,对求解得到的数值解进行可视化,以便对物质传输和扩散过程进行直观的理解。
三、MATLAB求解一维对流扩散方程的举例为了进一步说明MATLAB求解一维对流扩散方程的方法,下面举一个简单的例子进行说明。
一类对流-扩散方程源项反问题的数值解法
一类对流-扩散方程源项反问题是指求解一类对流-
扩散方程的源项,即求解源项函数$f(x,t)$,使得方程
$$\frac{\partial u}{\partial t}+\frac{\partial}{\partial x}(u\cdot
f(x,t))=0$$
的解满足给定的初始条件和边界条件。
解决一类对流-
扩散方程源项反问题的数值解法主要有以下几种:
(1)有限差分法:有限差分法是一种基于差分格式的数值解法,它将微分方程转化为一组线性方程组,然后使用数值求解方法求解。
(2)有限元法:有限元法是一种基于有限元的数值解法,它将微分方程转化为一组线性方程组,然后使用数值求解方法求解。
(3)有限体积法:有限体积法是一种基于有限体积的数值解法,它将微分方程转化为一组线性方程组,然后使用数值求解方法求解。
(4)有限元素法:有限元素法是一种基于有限元素的数值解法,它将微分方程转化为一组线性方程组,然后使用数值求解方法求解。
(5)积分变换法:积分变换法是一种基于积分变换的数值解法,它将微分方程转化为一组线性方程组,然后使用数值求解方法求解。
一维对流扩散方程的数值解法
对流-扩散方程是守恒定律控制方程的一种模型方程,它既是能量方程的表示形式,同时也可以认为是把压力梯度项隐含到了源项中去的动量方程的代表。
因此,以对流-扩散方程为例,来研究数值求解偏微分方程的相容性、收敛性和稳定性具有代表性的意义。
1 数学模型
本作业从最简单的模型方程,即一维、稳态、无源项的对流扩散方程出发,方程如下: 22, 02f f f
U D x t x x
∂∂∂+=≤≤∂∂∂ (1)
初始条件 (),0sin(2)f x t A kx π==
(2)
解析解
()()()224,sin 2Dk t
f x t e
A k x Ut ππ-=-
(3)
式中,1,0.05,0.5,1U D A k ====
函数(3)描述的是一个衰减波的图像,如图1所示
t=0 t=0.5 t=1
图1 函数()()()224,sin 2Dk t
f x t e
k x Ut ππ-=- 的图像(U=1,D=0.05,k=1)
2 数值解法
2.1 数值误差分析
在网格点(),i n 上差分方程的数值解n
i f 偏离该点上相应的偏微分方程的精确解
(),f i n 的值,称为网格节点上的数值误差。
当取定网格节点数21N =时,观察差分方程的解与微分方程的解在不同时间步长下的趋近程度,其中时间步长分别取值0.05,0.025,0.0125,0.0005t ∆=。
(a )21,0.05N t =∆= (b )21,0.025N t =∆=
(c )21,0.0125N t =∆= (d )201,0.0005N t =∆=
图2 数值误差随步长的变化情况
从图2的(a)~(d)可以定性的看出,数值误差与步长的大小有关。
在满足稳定性条件的前提下,数值误差随着时间步长的减小而减小,同时,图(d )表示增大网格的分辨率也有助于减小网格误差。
为了对数值误差有一个定量的认识,接下来取定时间步长为0.0005t ∆=,分别算出
11,21,41,61,81,101,121,161N =时,指标E =1所示。
表1 不同网格节点数下指标E 的值
由表1可以看出,指标E 的值随着网格节点数的增加呈先大幅减小后平缓变化的趋势,该趋势示于图3中。
(a ) (b )
图3 不同网格节点数下指标E 的值
在图3的(a )和(b )中,随着网格节点数的增加,曲线最后都趋于水平,即此时再增加网格节点数,对于提高数值求解的精确性作用不大。
2.2 截断误差分析
对于方程
()5f x x =
(4)
其一阶导数和二阶导数的精确解分别为4
5f x x
∂=∂和23220f x x ∂=∂。
一阶导数和二阶导数的中心差分表示式分别为
211()2i i f f f
o h x h
+--∂=+∂和221122
2()i i i f f f f
o h x h
+--+∂=+∂,2()o h 为截断误差。
当空间步长分别取值为0.2,0.1,0.05,0.025,0.0125h =时,节点()1i x =取处一阶导数和二阶导数的精确解与离散解的差值见表2。
表2 截断误差随空间步长的变化
图示如下
图4 截断误差随空间步长的变化
由图4可以看出,截断误差随着步长的减小而减小,当步长趋近于0时,截断误差趋近于0。
3 结论
(1)在满足稳定性条件的前提下,数值误差随着时间步长的减小而减小。
而对于空间步长的减小,数值误差呈先大幅减小后变平缓的趋势,即存在一个网格无关解的网格节点数。
(2)对于导数在空间上的离散,二阶精度的中心差分格式能较好的满足数值求解精确性的要求,当选取合适的步长时,其精确性也将显著提高。
%% =============================================================== clear
clc
format long
%% ===============================================================
U=1;D=0.05;k=1;A=0.5;
T=0.5;dt=0.05;m=T/dt;
L=2;N=21;h=L/(N-1);
%% ===============================================================
E0=0; E=0; %误差估计
%% ===============================================================
x=linspace(0,L,N);y=[1:N];
for i=1:N
f0(i)=A*sin(2*pi*k*x(i)); %初始条件
end
%% ===============================================================
for is=1:m,is
%% =============================================================== for i=2:N-1
f(i)=f0(i)-(U*dt/(2*h))*(f0(i+1)-f0(i-1))+(D*dt/(h^2))*(f0(i+1)-2 *f0(i)+f0(i-1));
end
%f(1)=f0(1)-(U*dt/(2*h))*(f0(2)-f0(N-1))+(D*dt/(h^2))*(f0(2)-2*f0 (1)+f0(N-1)); %边界条件
%f(N)=f0(N)-(U*dt/(2*h))*(f0(2)-f0(N-1))+(D*dt/(h^2))*(f0(2)-2*f0 (N)+f0(N-1)); %边界条件
f(1)=exp(-4*D*k^2*pi^2*dt*is)*A*sin(2*pi*k*(0-U*dt*is)); %边界条件 f(N)=exp(-4*D*k^2*pi^2*dt*is)*A*sin(2*pi*k*(L-U*dt*is)); %±边界条件 f0=f;
end
%% ===============================================================
for i=1:N
fexact(i)=exp(-4*D*k^2*pi^2*T)*A*sin(2*pi*k*(x(i)-U*T));
E0=(f(i)-fexact(i))^2+E0;
end
E=h*sqrt(E0)
%% ===============================================================
f1=exp(-4*D*k^2*pi^2*T)*A*sin(2*pi*k*(x-U*T));
plot(y,f1,'r',y,f,'b','linewidt',2),legend('Exact','Numerical'),axis( [1 N -2 2]),grid on
%plot(y,f1,'r','linewidt',4),legend('Exact'),axis([1 N -2 2]),grid on %hold on
%plot(y,f,'b','linewidt',2),legend('Numerical')
%==================================================================== i=1; %Çói(x=1)½ÚµãÉϵĵ¼ÊýÖµ
h=[0.2,0.1,0.05,0.025,0.0125];
%==================================================================== for n=1:5
f1(n)=((i+h(n))^5-(i-h(n))^5)/(2*h(n)); %f=x^5;
f2(n)=((i+h(n))^5-2*i^5+(i-h(n))^5)/(h(n))^2;
E1(n)=f1(n)-5*i^4;
E2(n)=f2(n)-20*i^2;
end
%==================================================================== loglog(1./h,E1,'o-',1./h,E2,'*-'),xlabel('1/h'),ylabel('Error'),... legend('first-order derivative','second-order derivative')。