有限差分法求解偏微分方程MATLAB教学教材

《使用 MATLAB 有限差分法求解非齐次偏微分方程》在科学和工程领域，偏微分方程是描述自然现象和过程中关键的数学工具。

非齐次偏微分方程作为其中的一个重要分支，在描述真实世界中的复杂现象方面具有广泛的应用。

而 MATLAB 作为一个强大的数学建模和计算工具，其有限差分法求解非齐次偏微分方程的能力受到了广泛关注。

在本文中，我们将以 MATLAB 为工具，探讨有限差分法如何用于求解非齐次偏微分方程，以及其中涉及的深度和广度。

1. 偏微分方程及有限差分法简介当我们研究自然界中的变化和现象时，经常会遇到连续变量之间的相关性和变化规律。

偏微分方程便是用来描述这些连续变量之间关系的数学工具。

而有限差分法则是一种数值计算方法，通过将连续的变量离散化，将偏微分方程转化为代数方程组，从而求解偏微分方程的数值解。

2. 非齐次偏微分方程的求解非齐次偏微分方程与常见的齐次偏微分方程相比，具有更复杂的边界和初始条件，因此其求解方法也更为复杂。

通过有限差分法，我们可以将非齐次偏微分方程转化为离散的代数方程组，进而求解出数值解。

3. MATLAB 中有限差分法的实现MATLAB 提供了丰富的数学建模和计算工具，包括用于求解偏微分方程的函数和工具箱。

通过调用这些函数和工具箱，我们可以方便地实现有限差分法对非齐次偏微分方程的求解。

4. 示例应用与个人观点我们将以一个实际的例子，展示 MATLAB 中有限差分法求解非齐次偏微分方程的过程，并共享对这一过程的个人观点和理解。

通过该示例，我们能更深刻地理解有限差分法在求解非齐次偏微分方程中的应用，以及其中涉及的数学原理和算法流程。

总结与回顾在本文中，我们以 MATLAB 为工具，探讨了有限差分法求解非齐次偏微分方程的深度和广度。

通过对有限差分法的基本原理和实际应用进行全面评估，我们详细介绍了有限差分法在求解非齐次偏微分方程中的具体步骤和流程。

我们也共享了在示例应用中对这一过程的个人理解和观点，以期帮助读者更全面、深刻和灵活地理解该主题。

matlab有限差分法求解非齐次偏微分方程【导语】本文将介绍matlab有限差分法在求解非齐次偏微分方程中的应用。

非齐次偏微分方程是数学和物理学中的常见问题之一，它们描述了许多实际系统的行为。

通过有限差分法，可以将偏微分方程转化为差分方程，从而利用计算机来求解。

本文将从原理、步骤和实例三个方面来分析非齐次偏微分方程的有限差分法求解过程。

【正文】一、原理有限差分法是将连续函数在一系列有限的点上进行逼近的方法。

它的基本思想是用差分代替微分，将偏导数转化为差分算子。

通过对空间和时间离散化，将非齐次偏微分方程转化为差分方程组，再利用数值计算的方法求解这个差分方程组，从而得到非齐次偏微分方程的近似解。

具体而言，有限差分法将求解区域划分为网格，并在网格上近似表示偏微分方程中的函数。

利用中心差分公式或向前、向后差分公式来近似计算偏导数。

通过将偏微分方程中的微分算子替换为差分近似，可以将方程转化为一个代数方程组，进而求解得到非齐次偏微分方程的近似解。

二、步骤1. 确定求解的区域和方程：首先要确定求解的区域，然后确定非齐次偏微分方程的形式。

在matlab中，可以通过定义一个矩阵来表示求解区域，并将方程转化为差分算子形式。

2. 离散化：将求解区域划分为网格，确定每个网格点的位置，建立网格点之间的连接关系。

通常，使用均匀网格来离散化求解区域，并定义网格点的坐标。

3. 建立差分方程组：根据偏微分方程的形式和离散化的结果，建立差分方程组。

根据中心差分公式，用网格点上的函数值和近邻点的函数值来近似计算偏导数。

将差分算子应用于非齐次偏微分方程的各个项，得到差分方程组。

4. 求解差分方程组：利用线性代数求解差分方程组。

将方程组转化为矩阵形式，利用matlab中的线性方程组求解功能，得到差分方程组的近似解。

通过调整求解区域划分的精细程度和差分算子的选取，可以提高求解的精度。

5. 回代和结果分析：将求解的结果回代到原非齐次偏微分方程中，分析其物理意义和数值稳定性。

南京理工大学课程考核论文课程名称：高等数值分析论文题目：有限差分法求解偏微分方程姓名：罗晨学号：成绩：有限差分法求解偏微分方程一、主要内容1.有限差分法求解偏微分方程，偏微分方程如一般形式的一维抛物线型方程：具体求解的偏微分方程如下：2.推导五种差分格式、截断误差并分析其稳定性；3.编写MATLAB程序实现五种差分格式对偏微分方程的求解及误差分析；4.结论及完成本次实验报告的感想。

二、推导几种差分格式的过程：有限差分法（finite-differencemethods ）是一种数值方法通过有限个微分方程近似求导从而寻求微分方程的近似解。

有限差分法的基本思想是把连续的定解区域用有限个离散点构成的网格来代替；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近似，积分用积分和来近似，于是原微分方程和定解条件就近似地代之以代数方程组，即有限差分方程组，解此方程组就可以得到原问题在离散点上的近似解。

推导差分方程的过程中需要用到的泰勒展开公式如下：()2100000000()()()()()()()......()(())1!2!!n n n f x f x f x f x f x x x x x x x o x x n +'''=+-+-++-+-(2-1)求解区域的网格划分步长参数如下：11k k k kt t x x h τ++-=⎧⎨-=⎩(2-2) 2.1古典显格式2.1.1古典显格式的推导由泰勒展开公式将(,)u x t 对时间展开得2,(,)(,)()()(())i i k i k k k uu x t u x t t t o t t t∂=+-+-∂(2-3) 当1k t t +=时有21,112,(,)(,)()()(())(,)()()i k i k i k k k k k i k i k uu x t u x t t t o t t tuu x t o tττ+++∂=+-+-∂∂=+⋅+∂(2-4)得到对时间的一阶偏导数1,(,)(,)()=()i k i k i k u x t u x t uo t ττ+-∂+∂(2-5) 由泰勒展开公式将(,)u x t 对位置展开得223,,21(,)(,)()()()()(())2!k i k i k i i k i i u uu x t u x t x x x x o x x x x∂∂=+-+-+-∂∂(2-6)当11i i x x x x +-==和时，代入式(2-6)得2231,1,1122231,1,1121(,)(,)()()()()(())2!1(,)(,)()()()()(())2!i k i k i k i i i k i i i i i k i k i k i i i k i i i iu u u x t u x t x x x x o x x x xu u u x t u x t x x x x o x x x x ++++----⎧∂∂=+-+-+-⎪⎪∂∂⎨∂∂⎪=+-+-+-⎪∂∂⎩(2-7) 因为1k k x x h +-=，代入上式得2231,,22231,,21(,)(,)()()()2!1(,)(,)()()()2!i k i k i k i k i k i k i k i ku u u x t u x t h h o h x xu u u x t u x t h h o h x x +-⎧∂∂=+⋅+⋅+⎪⎪∂∂⎨∂∂⎪=-⋅+⋅+⎪∂∂⎩(2-8) 得到对位置的二阶偏导数2211,22(,)2(,)(,)()()i k i k i k i k u x t u x t u x t uo h x h+--+∂=+∂(2-9) 将式(2-5)、(2-9)代入一般形式的抛物线型偏微分方程得(2-10)为了方便我们可以将式(2-10)写成11122k kk k k k i i i i i i u u u u u f h ατ++-⎡⎤--+-=⎢⎥⎣⎦(2-11) ()11122k k k k k k i i i i i i u u uu u f hτατ++----+=(2-12)最后得到古典显格式的差分格式为()111(12)k k k k k i i i i i u ra u r u u f ατ++-=-+++(2-13)2r hτ=其中，古典显格式的差分格式的截断误差是2()o h τ+。

引言偏微分方程定解问题有着广泛的应用背景。

人们用偏微分方程来描述、解释或者预见各种自然现象，并用于科学和工程技术的各个领域fll。

然而，对于广大应用工作者来说，从偏微分方程模型出发，使用有限元法或有限差分法求解都要耗费很大的工作量，才能得到数值解。

现在，MATLAB PDEToolbox已实现对于空间二维问题高速、准确的求解过程。

偏微分方程如果一个微分方程中出现的未知函数只含一个自变量，这个方程叫做常微分方程，也简称微分方程；如果一个微分方程中出现多元函数的偏导数，或者说如果未知函数和几个变量有关，而且方程中出现未知函数对几个变量的导数，那么这种微分方程就是偏微分方程。

常用的方法有变分法和有限差分法。

变分法是把定解问题转化成变分问题，再求变分问题的近似解；有限差分法是把定解问题转化成代数方程，然后用计算机进行计算；还有一种更有意义的模拟法，它用另一个物理的问题实验研究来代替所研究某个物理问题的定解。

虽然物理现象本质不同，但是抽象地表示在数学上是同一个定解问题，如研究某个不规则形状的物体里的稳定温度分布问题，由于求解比较困难，可作相应的静电场或稳恒电流场实验研究，测定场中各处的电势，从而也解决了所研究的稳定温度场中的温度分布问题。

随着物理科学所研究的现象在广度和深度两方面的扩展，偏微分方程的应用范围更广泛。

从数学自身的角度看，偏微分方程的求解促使数学在函数论、变分法、级数展开、常微分方程、代数、微分几何等各方面进行发展。

从这个角度说，偏微分方程变成了数学的中心。

一、MATLAB方法简介及应用1.1 MATLAB简介MATLAB是美国MathWorks公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MATLAB和Simulink两大部分。

1.2 Matlab主要功能数值分析数值和符号计算工程与科学绘图控制系统的设计与仿真数字图像处理数字信号处理通讯系统设计与仿真财务与金融工程1.3 优势特点1) 高效的数值计算及符号计算功能，能使用户从繁杂的数学运算分析中解脱出来；2) 具有完备的图形处理功能，实现计算结果和编程的可视化；3) 友好的用户界面及接近数学表达式的自然化语言，使学者易于学习和掌握；4) 功能丰富的应用工具箱(如信号处理工具箱、通信工具箱等) ，为用户提供了大量方便实用的处理工具。

引言偏微分方程定解问题有着广泛的应用背景。

人们用偏微分方程来描述、解释或者预见各种自然现象，并用于科学和工程技术的各个领域fll。

然而，对于广大应用工作者来说，从偏微分方程模型出发，使用有限元法或有限差分法求解都要耗费很大的工作量，才能得到数值解。

现在，MATLAB PDEToolbox已实现对于空间二维问题高速、准确的求解过程。

偏微分方程如果一个微分方程中出现的未知函数只含一个自变量，这个方程叫做常微分方程，也简称微分方程；如果一个微分方程中出现多元函数的偏导数，或者说如果未知函数和几个变量有关，而且方程中出现未知函数对几个变量的导数，那么这种微分方程就是偏微分方程。

常用的方法有变分法和有限差分法。

变分法是把定解问题转化成变分问题，再求变分问题的近似解；有限差分法是把定解问题转化成代数方程，然后用计算机进行计算；还有一种更有意义的模拟法，它用另一个物理的问题实验研究来代替所研究某个物理问题的定解。

虽然物理现象本质不同，但是抽象地表示在数学上是同一个定解问题，如研究某个不规则形状的物体里的稳定温度分布问题，由于求解比较困难，可作相应的静电场或稳恒电流场实验研究，测定场中各处的电势，从而也解决了所研究的稳定温度场中的温度分布问题。

随着物理科学所研究的现象在广度和深度两方面的扩展，偏微分方程的应用范围更广泛。

从数学自身的角度看，偏微分方程的求解促使数学在函数论、变分法、级数展开、常微分方程、代数、微分几何等各方面进行发展。

从这个角度说，偏微分方程变成了数学的中心。

一、MATLAB方法简介及应用1.1 MATLAB简介MATLAB是美国MathWorks公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MATLAB和Simulink两大部分。

1.2 Matlab主要功能数值分析数值和符号计算工程与科学绘图控制系统的设计与仿真数字图像处理数字信号处理通讯系统设计与仿真财务与金融工程1.3 优势特点1) 高效的数值计算及符号计算功能，能使用户从繁杂的数学运算分析中解脱出来；2) 具有完备的图形处理功能，实现计算结果和编程的可视化；3) 友好的用户界面及接近数学表达式的自然化语言，使学者易于学习和掌握；4) 功能丰富的应用工具箱(如信号处理工具箱、通信工具箱等) ，为用户提供了大量方便实用的处理工具。

南京理工大学课程考核论文课程名称：高等数值分析论文题目：有限差分法求解偏微分方程*名：**学号： 1成绩：有限差分法求解偏微分方程一、主要内容1.有限差分法求解偏微分方程，偏微分方程如一般形式的一维抛物线型方程：22(,)()u uf x t t xαα∂∂-=∂∂其中为常数具体求解的偏微分方程如下：22001(,0)sin()(0,)(1,)00u u x t x u x x u t u t t π⎧∂∂-=≤≤⎪∂∂⎪⎪⎪=⎨⎪⎪==≥⎪⎪⎩2.推导五种差分格式、截断误差并分析其稳定性；3.编写MATLAB 程序实现五种差分格式对偏微分方程的求解及误差分析；4.结论及完成本次实验报告的感想。

二、推导几种差分格式的过程：有限差分法（finite-difference methods ）是一种数值方法通过有限个微分方程近似求导从而寻求微分方程的近似解。

有限差分法的基本思想是把连续的定解区域用有限个离散点构成的网格来代替；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近似，积分用积分和来近似，于是原微分方程和定解条件就近似地代之以代数方程组，即有限差分方程组，解此方程组就可以得到原问题在离散点上的近似解。

推导差分方程的过程中需要用到的泰勒展开公式如下：()2100000000()()()()()()()......()(())1!2!!n n n f x f x f x f x f x x x x x x x o x x n +'''=+-+-++-+- (2-1)求解区域的网格划分步长参数如下：11k k k kt t x x h τ++-=⎧⎨-=⎩ (2-2) 2.1 古典显格式2.1.1 古典显格式的推导由泰勒展开公式将(,)u x t 对时间展开得 2,(,)(,)()()(())i i k i k k k uu x t u x t t t o t t t∂=+-+-∂ (2-3) 当1k t t +=时有21,112,(,)(,)()()(())(,)()()i k i k i k k k k k i k i k uu x t u x t t t o t t tuu x t o tττ+++∂=+-+-∂∂=+⋅+∂ (2-4)得到对时间的一阶偏导数1,(,)(,)()=()i k i k i k u x t u x t uo t ττ+-∂+∂ (2-5) 由泰勒展开公式将(,)u x t 对位置展开得223,,21(,)(,)()()()()(())2!k i k i k i i k i i u uu x t u x t x x x x o x x x x∂∂=+-+-+-∂∂ (2-6)当11i i x x x x +-==和时，代入式(2-6)得2231,1,1122231,1,1121(,)(,)()()()()(())2!1(,)(,)()()()()(())2!i k i k i k i i i k i i i i i k i k i k i i i k i i i iu uu x t u x t x x x x o x x x x u u u x t u x t x x x x o x x x x ++++----⎧∂∂=+-+-+-⎪⎪∂∂⎨∂∂⎪=+-+-+-⎪∂∂⎩(2-7) 因为1k k x x h +-=，代入上式得2231,,22231,,21(,)(,)()()()2!1(,)(,)()()()2!i k i k i k i k i k i k i k i ku uu x t u x t h h o h x xu u u x t u x t h h o h x x +-⎧∂∂=+⋅+⋅+⎪⎪∂∂⎨∂∂⎪=-⋅+⋅+⎪∂∂⎩ (2-8) 得到对位置的二阶偏导数2211,22(,)2(,)(,)()()i k i k i k i k u x t u x t u x t u o h x h+--+∂=+∂ (2-9) 将式(2-5)、(2-9)代入一般形式的抛物线型偏微分方程得21112(,)(,)(,)2(,)(,)(,)()i k i k i k i k i k i k u x t u x t u x t u x t u x t f x t o h h αττ++---+⎡⎤-=++⎢⎥⎣⎦(2-10)为了方便我们可以将式(2-10)写成11122k kk k k k i i i i i i u u u u u f h ατ++-⎡⎤--+-=⎢⎥⎣⎦(2-11) ()11122k k k k k k i i i i i i u u uu u f h τατ++----+= (2-12)最后得到古典显格式的差分格式为()111(12)k k k k k i i i i i u ra u r u u f ατ++-=-+++ (2-13)2r h τ=其中，古典显格式的差分格式的截断误差是2()o h τ+。

《偏微分方程数值解》上机报告实验内容 1:分别用向前差分格式、向后差分格式及六点对称格式, 求解下列问题:222, 01, 0, (0, (1, 0, 1, (, 0 sin( (1.u u x t t x u t u t t u x x x x π⎧∂∂=+<<>⎪∂∂⎪⎨==>⎪⎪=+−⎩x 方向 0.1h =, t 方向0.01τ=.在 0.25t =时观察数值解与精确解 2sin( (1 u e x x x ππ−=+−的误差. (一算法描述:(二实验结果:1.误差的数值解结果数值对比(A“向前差分格式”程序:>>forward(0.1,0.01, 0.25Current plot heldans =0.00000.00270.00510.00700.00820.00870.00820.00700.00510.00270.0000(B“向后差分格式”程序:>>back(0.1,0.01, 0.25Current plot heldans =0.0000-0.0037-0.0071-0.0097-0.0114-0.0120-0.0114-0.0097-0.0071 -0.00370.0000(C“六点差分格式”程序:>>six(0.1,0.01, 0.25Current plot heldans =0.0000-0.0005-0.0009-0.0013-0.0015-0.0016-0.0015-0.0013-0.0009-0.00050.0000注:这里的"误差"=精确解-数值解.2.精确解与数值解结果图像对比“向前差分格式”:注:曲线表示精确解,"o"表示数值解(t=0.25时. “向后差分格式”:注:曲线表示精确解,"o"表示数值解(t=0.25时. “六点差分格式” :注:曲线表示精确解,"O"表示数值解(t=0.25时.(三结果分析通过(一 , (二 ,我们检验了三种方法都能很好的求解此一维热传导方程,其中明显能发现“六点对称格式”的误差更小。

第四讲【2 】Matlab求解微分方程（组）理论介绍：Matlab求解微分方程（组）敕令求解实例：Matlab求解微分方程（组）实例现实运用问题经由过程数学建模所归纳得到的方程,绝大多半都是微分方程,真正能得到代数方程的机遇很少.另一方面,可以或许求解的微分方程也是十分有限的,特别是高阶方程和偏微分方程（组）.这就请求我们必须研讨微分方程（组）的解法：解析解法和数值解法.一．相干函数.敕令及简介1.在Matlab中,用大写字母D表示导数,Dy表示y关于自变量的一阶导数,D2y表示y关于自变量的二阶导数,依此类推.函数dsolve用来解决常微分方程（组）的求解问题,挪用格局为：X=dsolve(‘eqn1’,’eqn2’,…)函数dsolve用来解符号常微分方程.方程组,假如没有初始前提,则求出通解,假如有初始前提,则求出特解.留意,体系缺省的自变量为t2.函数dsolve求解的是常微分方程的准确解法,也称为常微分方程的符号解.但是,有大量的常微分方程固然从理论上讲,其解是消失的,但我们却无法求出其解析解,此时,我们须要追求方程的数值解,在求常微分方程数值解方面,MATLAB 具有丰硕的函数,我们将其统称为solver,其一般格局为：[T,Y]=solver(odefun,tspan,y0)解释：(1)solver为敕令ode45.ode23.ode113.ode15s.ode23s.ode23t.ode23tb.ode15i之一.(2)odefun 是显示微分方程'(,)y f t y =在积分区间tspan 0[,]f t t =上从0t 到f t 用初始前提y 求解.(3)假如要获得微分方程问题在其他指准时光点012,,,,ft t t t 上的解,则令tspan012[,,,]f t t t t =(请求是单调的).(4)因为没有一种算法可以有用的解决所有的ODE 问题,为此,Matlab 供给了多种求解器solver,对于不同的ODE 问题,采用不同的solver.表1 Matlab 中文本文件读写函数解释：ode23.ode45是极其常用的用来求解非刚性的标准情势的一阶微分方程（组）的初值问题的解的Matlab 常用程序,个中：ode23采用龙格-库塔2阶算法,用3阶公式作误差估量来调节步长,具有低等的精度.ode45则采用龙格-库塔4阶算法,用5阶公式作误差估量来调节步长,具有中等的精度.3．在matlab敕令窗口.程序或函数中创建局部函数时,可用内联函数inline,inli ne函数情势相当于编写M函数文件,但不需编写M-文件就可以描写出某种数学关系.挪用inline函数,只能由一个matlab表达式构成,并且只能返回一个变量,不许可[u,v]这种向量情势.因而,任何请求逻辑运算或乘法运算以求得最终成果的场合,都不能运用inline函数,inline函数的一般情势为：FunctionName=inline(‘函数内容’, ‘所有自变量列表’)例如：（求解F(x)=x^2*cos(a*x)-b ,a,b是标量;x是向量）在敕令窗口输入: Fofx=inline(‘x .^2*cos(a*x)-b’ , ‘x’,’a’,’b’);g= Fofx([pi/3 pi/3.5],4,1)体系输出为：g=-1.5483 -1.7259留意：因为运用内联对象函数inline不须要别的树立m文件,所有运用比较便利,别的在运用ode45函数的时刻,界说函数往往须要编辑一个m文件来单独界说,如许不便于治理文件,这里可以运用inline来界说函数.二．实例介绍1.几个可以直接用Matlab求微分方程准确解的实例例1 求解微分方程2 '2x y xy xe-+=程序：syms x y; y=dsolve(‘Dy+2*x*y=x*exp(-x^2)’,’x’)例 2求微分方程'0xxy y e+-=在初始前提(1)2y e=下的特解并画出解函数的图形.程序：syms x y; y=dsolve(‘x*Dy+y-exp(1)=0’,’y(1)=2*exp(1)’,’x’);ezplot(y)例 3求解微分方程组530tdx x y e dt dy x y dt ⎧++=⎪⎪⎨⎪--=⎪⎩在初始前提00|1,|0t t x y ====下的特解并画出解函数的图形.程序：syms x y t[x,y]=dsolve('Dx+5*x+y=exp(t)','Dy-x-3*y=0','x(0)=1','y(0)=0','t') simple(x); simple(y)ezplot(x,y,[0,1.3]);axis auto2.用ode23.ode45等求解非刚性标准情势的一阶微分方程（组）的初值问题的数值解（近似解）例4求解微分方程初值问题2222(0)1dy y x x dx y ⎧=-++⎪⎨⎪=⎩的数值解,求解规模为区间[0,0.5].程序：fun=inline('-2*y+2*x^2+2*x','x','y'); [x,y]=ode23(fun,[0,0.5],1); plot(x,y,'o-') 例5求解微分方程22'2(1)0,(0)1,(0)0d y dy y y y y dt dt μ--+===的解,并画出解的图形.剖析：这是一个二阶非线性方程,我们可以经由过程变换,将二阶方程化为一阶方程组求解.令12,,7dyx y x dt μ===,则121221212,(0)17(1),(0)0dx x x dtdx x x x x dt ⎧==⎪⎪⎨⎪=--=⎪⎩编写M-文件vdp.m function fy=vdp(t,x)fy=[x(2);7*(1-x(1)^2)*x(2)-x(1)]; end在Matlab 敕令窗口编写程序 y0=[1;0][t,x]=ode45(@vdp,[0,40],y0);或[t,x]=ode45('vdp',[0,40],y0); y=x(:,1);dy=x(:,2); plot(t,y,t,dy)演习与思虑：M-文件vdp.m 改写成inline 函数程序？ 3.用Euler 折线法求解Euler 折线法求解的根本思惟是将微分方程初值问题00(,)()dyf x y dx y x y ⎧=⎪⎨⎪=⎩化成一个代数(差分)方程,重要步骤是用差商()()y x h y x h +-替代微商dydx ,于是00()()(,())()k k k k y x h y x f x y x h y y x +-⎧=⎪⎨⎪=⎩记1,(),k k k k x x h y y x +=+=从而1(),k k y y x h +=+于是0011(),,0,1,2,,1(,).k k k k k k y y x x x h k n y y hf x y ++=⎧⎪=+=-⎨⎪=+⎩例6用Euler 折线法求解微分方程初值问题22(0)1dyx y dxy y ⎧=+⎪⎨⎪=⎩的数值解（步长h 取0.4）,求解规模为区间[0,2].剖析：本问题的差分方程为00110,1,0.4,0,1,2,,1(,).k k k k k k x y h x x h k n y y hf x y ++===⎧⎪=+=-⎨⎪=+⎩程序：>> clear >> f=sym('y+2*x/y^2'); >> a=0; >> b=2; >> h=0.4; >> n=(b-a)/h+1; >> x=0; >> y=1;>> szj=[x,y];%数值解 >> for i=1:n-1y=y+h*subs(f,{'x','y'},{x,y});%subs,调换函数 x=x+h; szj=[szj;x,y];end >>szj>> plot(szj(:,1),szj(:,2)) 解释：调换函数subs 例如：输入subs(a+b,a,4)意思就是把a 用4调换掉落,返回4+b,也可以调换多个变量,例如：subs(cos(a)+sin(b),{a,b},[sym('alpha'),2])分离用字符alpha 调换a 和2调换b,返回 cos(alpha)+sin(2)特别解释：本问题可进一步运用四阶Runge-Kutta 法求解,Euler 折线法现实上就是一阶Runge-Kutta 法,Runge-Kutta 法的迭代公式为001112341213243(),,(22),6(,),0,1,2,,1(,),22(,),22(,).k k k k k k k k k k k k y y x x x h h y y L L L L L f x y k n h h L f x y L h h L f x y L L f x h y hL ++=⎧⎪=+⎪⎪=++++⎪⎪=⎪=-⎨⎪=++⎪⎪⎪=++⎪⎪=++⎩响应的Matlab 程序为：>> clear >> f=sym('y+2*x/y^2'); >> a=0; >> b=2; >> h=0.4; >> n=(b-a)/h+1;>> x=0; >> y=1;>> szj=[x,y];%数值解 >> for i=1:n-1l1=subs(f,{'x','y'},{x,y});调换函数 l2=subs(f,{'x','y'},{x+h/2,y+l1*h/2}); l3=subs(f,{'x','y'},{x+h/2,y+l2*h/2}); l4=subs(f,{'x','y'},{x+h,y+l3*h}); y=y+h*(l1+2*l2+2*l3+l4)/6; x=x+h; szj=[szj;x,y]; end >>szj>> plot(szj(:,1),szj(:,2))演习与思虑：(1)ode45求解问题并比较差异.(2)运用Matlab 求微分方程(4)(3)''20y y y -+=的解. (3)求解微分方程''2',2(1)0,030,(0)1,(0)0y y y y x y y --+=≤≤==的特解.(4)运用Matlab 求微分方程初值问题2''''00(1)2,|1,|3x x x y xy y y ==+===的解.提示：尽可能多的斟酌解法三．微分方程转换为一阶显式微分方程组Matlab 微分方程解算器只能求解标准情势的一阶显式微分方程（组）问题,是以在运用ODE 解算器之前,我们须要做的第一步,也是最重要的一步就是借助状况变量将微分方程（组）化成Matlab 可接收的标准情势.当然,假如ODEs 由一个或多个高阶微分方程给出,则我们应先将它变换成一阶显式常微分方程组.下面我们以两个高阶微分方程组构成的ODEs 为例介绍若何将它变换成一个一阶显式微分方程组.Step1将微分方程的最高阶变量移到等式左边,其它移到右边,并按阶次从低到高分列.情势为：()'''(1)'''(1)()'''(1)'''(1)(,,,,,,,,,,)(,,,,,,,,,,)m m n n m n x f t x x x x y y y y y g t x x x xy y y y ----⎧=⎨=⎩ Step 2 为每一阶微分式选择状况变量,最高阶除外'''(1)123'''(1)123,,,,,,,,,m m n m m m m n x x x x x x x x x y x y x y x y --++++========留意：ODEs 中所有是因变量的最高阶次之和就是须要的状况变量的个数,最高阶的微分式不须要给它状况变量.Step 3 依据选用的状况变量,写出所有状况变量的一阶微分表达式''''122334123''12123,,,,(,,,,,),,(,,,,,)m m n m m m n m n x x x x x x x f t x x x x x x x g t x x x x +++++======演习与思虑：(1)求解微分方程组**'''3312*'''3312()()22x x x y x r r y y y x y r r μμμμμμ⎧+-=+--⎪⎪⎨⎪=+--⎪⎩个中2r =1r =*1,μμ=-1/82.45,μ=(0) 1.2,x =(0)0,y ='(0)0,x ='(0) 1.049355751y =-(2)求解隐式微分方程组''''''''''''2235x y x y x y x y xy y ⎧+=⎨++-=⎩提示：运用符号盘算函数solve 求'''',x y ,然后运用求解微分方程的办法四．偏微分方程解法Matlab 供给了两种办法解决PDE 问题,一是运用pdepe 函数,它可以求解一般的PDEs,具有较大的通用性,但只支撑敕令情势挪用;二是运用PDE 对象箱,可以求解特别PDE 问题,PDEtoll 有较大的局限性,比如只能求解二阶PDE 问题,并且不能解决片微分方程组,但是它供给了GUI 界面,从庞杂的编程中摆脱出来,同时还可以经由过程File —>Save As 直接生成M 代码.1.一般偏微分方程（组）的求解(1)Matlab 供给的pdepe 函数,可以直接求解一般偏微分方程（组）,它的挪用格局为：sol=pdepe(m,@pdefun,@pdeic,@pdebc,x,t)@pdefun 是PDE 的问题描写函数,它必须换成标准情势：(,,)[(,,,)](,,,)m m u u u uc x t x x f x t u s x t u x t x x x -∂∂∂∂∂=+∂∂∂∂∂如许,PDE 就可以编写进口函数：[c,f,s]=pdefun(x,t,u,du),m,x,t 对应于式中相干参数,du 是u 的一阶导数,由给定的输入变量可表示出c,f,s 这三个函数.@pdebc 是PDE 的边界前提描写函数,它必须化为情势：(,,)(,,).*(,,,)0up x t u q x t u f x t u x ∂==∂于是边值前提可以编写函数描写为：[pa,qa,pb,qb]=pdebc(x,t,u,du),个中a 表示下边界,b 表示上边界.@pdeic 是PDE 的初值前提,必须化为情势：00(,)u x t u =,故可以运用函数描写为：u0=pdeic(x)sol 是一个三维数组,sol(:,:,i)表示i u 的解,换句话说,k u 对应x(i)和t(j)时的解为s ol(i,j,k),经由过程sol,我们可以运用pdeval 函数直接盘算某个点的函数值.(2)实例解释求解偏微分2111222221220.024()0.17()u u F u u t x u u F u u t x ⎧∂∂=--⎪⎪∂∂⎨∂∂⎪=+-⎪∂∂⎩个中, 5.7311.46()x x F x e e -=-且知足初始前提12(,0)1,(,0)0u x u x ==及边界前提1(0,)0,u t x ∂=∂221(0,)0,(1,)1,(1,)0u u t u t t x ∂===∂解：(1)对比给出的偏微分方程和pdepe 函数求解的标准情势,原方程改写为111221220.024()1.*()10.17u u F u u x u F u u u t x x ∂⎡⎤⎢⎥--⎡⎤⎡⎤⎡⎤∂∂∂=+⎢⎥⎢⎥⎢⎥⎢⎥-∂∂∂⎣⎦⎣⎦⎣⎦⎢⎥⎢⎥∂⎣⎦ 可见1121220.024()10,,,()10.17u F u u x m c f s F u u u x ∂⎡⎤⎢⎥--⎡⎤⎡⎤∂====⎢⎥⎢⎥⎢⎥-∂⎣⎦⎣⎦⎢⎥⎢⎥∂⎣⎦%目的PDE 函数function [c,f,s]=pdefun(x,t,u,du)c=[1;1];f=[0.024*du(1);0.17*du(2)];temp=u(1)-u(2);s=[-1;1].*(exp(5.73*temp)-exp(-11.46*temp))end(2)边界前提改写为：下边界2010.*00f u ⎡⎤⎡⎤⎡⎤+=⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦上边界1110.*000u f -⎡⎤⎡⎤⎡⎤+=⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦ %边界前提函数function [pa,qa,pb,qb]=pdebc(xa,ua,xb,ub,t)pa=[0;ua(2)];qa=[1;0];pb=[ub(1)-1;0];qb=[0;1];end(3)初值前提改写为：1210u u ⎡⎤⎡⎤=⎢⎥⎢⎥⎣⎦⎣⎦%初值前提函数function u0=pdeic(x)u0=[1;0];end(4)编写主调函数clcx=0:0.05:1;t=0:0.05:2;m=0;sol=pdepe(m,@pdefun,@pdeic,@pdebc,x,t);subplot(2,1,1)surf(x,t,sol(:,:,1))subplot(2,1,2)surf(x,t,sol(:,:,2))演习与思虑： This example illustrates the straightforward formulation, computation, and plotting of the solution of a single PDE.2()u u t x x π∂∂∂=∂∂∂This equation holds on an interval 01x ≤≤ for times 0t ≥. The PDE satisfies the initial condition (,0)sin u x x π= and boundary conditions(0,)0;(1,)0t u u t e t x π-∂=+=∂2.PDEtool 求解偏微分方程(1)PDEtool （GUI ）求解偏微分方程的一般步骤在Matlab 敕令窗口输入pdetool,回车,PDE 对象箱的图形用户界面(GUI)体系就启动了.从界说一个偏微分方程问题到完成解偏微分方程的定解,全部进程大致可以分为六个阶段Step 1 “Draw 模式”绘制平面有界区域Ω,经由过程公式把Matlab 体系供给的实体模子：矩形.圆.椭圆和多边形,组合起来,生成须要的平面区域.Step 2 “Boundary 模式”界说边界,声明不同边界段的边界前提.Step3 “PDE 模式”界说偏微分方程,肯定方程类型和方程系数c,a,f,d,依据具体情况,还可以在不同子区域声明不同系数.Step 4 “Mesh 模式”网格化区域Ω,可以掌握主动生成网格的参数,对生成的网格进行多次细化,使网格朋分更细更合理.Step 5 “Solve 模式”解偏微分方程,对于椭圆型方程可以激活并掌握非线性自顺应解题器来处理非线性方程;对于抛物线型方程和双曲型方程,设置初始边界前提后可以求出给准时刻t 的解;对于特点值问题,可以求出给定区间上的特点值.求解完成后,可以返回到Step 4,对网格进一步细化,进行再次求解.Step 6 “View 模式”盘算成果的可视化,可以经由过程设置体系供给的对话框,显示所求的解的表面图.网格图.等高线图和箭头梯形图.对于抛物线型和双曲线型问题的解还可以进行为画演示.(2)实例解释用法求解一个正方形区域上的特点值问题：12|0u u u u λ∂Ω⎧-∆-=⎪⎨⎪=⎩正方形区域为：11,1 1.x x -≤≤-≤≤(1)运用PDE 对象箱打开GUI 求解方程(2)进入Draw 模式,绘制一个矩形,然后双击矩形,在弹出的对话框中设置Left=-1,Bottom=-1,Width=2,Height=2,确认并封闭对话框(3)进入Boundary 模式,边界前提采用Dirichlet 前提的默认值(4)进入PDE 模式,单击对象栏PDE 按钮,在弹出的对话框中方程类型选择Eigenmodes,参数设置c=1,a=-1/2,d=1,确认后封闭对话框(5)单击对象栏的 按钮,对正方形区域进行初始网格剖分,然后再对网格进一步细化剖分一次(6)点开solve菜单,单击Parameters选项,在弹出的对话框中设置特点值区域为[ -20,20](7)单击Plot菜单的Parameters项,在弹出的对话框中选中Color.Height(3-D plot)和show mesh项,然后单击Done确认(8)单击对象栏的“=”按钮,开端求解。