Matlab解微分方程

第十六章 偏微分方程的数值解法

科学研究和工程技术中的许多问题可建立偏微分方程的数学模型。包含多个自变量的微分方程称为偏微分方程(partial differential equation)，简称PDE 。偏微分方程问题，其求解是十分困难的。除少数特殊情况外，绝大多数情况均难以求出精确解。因此，近似解法就显得更为重要。本章仅介绍求解各类典型偏微分方程定解问题的差分方法。

16.1 几类偏微分方程的定解问题

一个偏微分方程的表示通常如下：

(,,,,)xx xy yy x y A B C f x y Φ+Φ+Φ=ΦΦΦ (16.1.1) 式中，,,A B C 是常数，称为拟线性(quasilinear)数。通常，存在3种拟线性方程： 双曲型(hyperbolic)方程：240B AC ->； 抛物线型(parabolic)方程：240B AC -=； 椭圆型(ellliptic)方程：240B AC -<。

16.1.2 双曲型方程

最简单形式为一阶双曲型方程：

0u u a t x

??+=?? (16.1.2) 物理中常见的一维振动与波动问题可用二阶波动方程：

22222u u a t x

??=?? (16.1.3) 描述，它是双曲型方程的典型形式。方程的初值问题为：

222220

0,(,0)()()t u u

a

t x t

x u x x u x x t ?ψ=???=>-∞<<+∞

?????

=??

??=-∞<<+∞

??? (16.1.4)

边界条件一般有三类，最简单的初边值问题为：

22222120

00,0(,0)(0,)(),(,)()0()t u u

a t T x l t x u x l

u t g t u l t g t t T u

x x t ?ψ=???==<<<

16.1.3 抛物型方程

其最简单的形式为一维热传导方程：

220(0)u u

a a t x

??-=>?? (16.1.8) 方程可以有两种不同类型的定解问题：

(1) 初值问题：

2200,(,0)()u u

a t x t x

u x x x ????-=>-∞<<+∞?

????=-∞<<+∞

?

(16.1.6)

(2) 初边值问题：

2212

00,0(,0)()0(0,)(),(,)()0u u

a t T x l t x u x x x l u t g t u l t g t t T

????-=<<<

=≤≤??==≤≤???

(16.1.7) 其中()x ?，1()g t ，2()g t 为已知函数，且满足连接条件：

12(0)(0),()(0)g l g ??== (16.1.8)

边界条件12(0,)(),(,)()u t g t u l t g t ==为第一类边界条件。

第二类和第三类边界条件为：

10

122()()

()()

x x l

u t u g t x u t u g t x λλ==???-=?????

??

?+=?????

0t T ≤≤ (16.1.9)

其中1()0t λ≥，2()0t λ≥。当12()()0t t λλ=≡时，为第二类边界条件，否则称为第三类边界条件。

16.1.4 椭圆型方程

其最典型、最简单的形式是泊松(Poisson)方程

2222(,)u u

u f x y x y

???=+=?? (16.1.10)

特别地，当(,)0f x y ≡时，即为拉普拉斯(Laplace)方程，又称为调和方程：

02222=??+??=?y

u

x u u (16.1.11)

Poisson 方程的第一边值问题为：

2222(,)(,)(,)(,)

(,)x y u u

f x y x y x y u x y x y ?∈Γ???+=∈Ω?????=Γ=?Ω

?

(16.1.12)

其中Ω为以Γ为边界的有界区域，Γ为分段光滑曲线，ΩΓ称为定解区域，(,)f x y ，(,)x y ?分别为Ω，Γ上的已知连续函数。

第二类和第三类边界条件可统一表示为：

(,)(,)x y u u x y α?∈Γ

??

?+= ?

???

n (16.1.13)

其中n 为边界Γ的外法线方向。当0α=时为第二类边界条件，0α≠时为第三类边界条件。

16.2 差分方法的基本概念

差分方法又称为有限差分方法或网格法，是求偏微分方程定解问题的数值解中应用最广泛的方法之一。

它的基本思想是：先对求解区域作网格剖分，将自变量的连续变化区域用有限离散点(网格点)集代替；将问题中出现的连续变量的函数用定义在网格点上离散变量的函数代替；通过用网格点上函数的差商代替导数，将含连续变量的偏微分方程定解问题化成只含有限个未知数的代数方程组(称为差分格式)。如果差分格式有解，且当网格无限变小时其解收敛于原微分方程定解问题的解，则差分格式的解就作为原问题的近似解(数值解)。

因此，用差分方法求偏微分方程定解问题一般需要解决以下问题： (1) 选取网格；

(2) 对微分方程及定解条件选择差分近似，列出差分格式； (3) 求解差分格式；

(4) 讨论差分格式解对于微分方程解的收敛性及误差估计。

下面，用一个简单的例子来说明用差分方法求解偏微分方程问题的一般过程及差分方法的基本概念。

设有一阶双曲型方程初值问题。

00,(,0)()u

u a t x t

x u x x ????+=>-∞<<+∞

?????=? (16.2.1)

选取网格：

图16.2.1 差分示意图

首先对定解区域{(,),0}D x t x t =-∞<<+∞≥作网格剖分，最简单常用一种网格是用两族分别平行于x 轴与t 轴的等距直线：

k x x kh ==, (0,1,2

,0,1,2,)j t t j k j τ===±±= (16.2.2)

将D 分成许多小矩形区域。这些直线称为网格线，其交点称为网格点，也称为节点，h 和τ分别称作x 方向和t 方向的步长。这种网格称为矩形网格。

(1) 对微分方程及定解条件选择差分近似，列出差分格式。如果用向前差商表示一阶偏导数，即：

2211(,)12(,)

(,)(,)(,)2(,)(,)(,)

2

k j k j k j k j k j x x t k j k j k j t x t u x t u x t u

h u x h t x h u x t u x t u u x t t θτθττ++-?'

'=-+?-?''=-+? (16.2.3)

其中120,1θθ<<。

方程：

0u u

a t x

??+=?? (16.2.4) 在节点(,)k j x t 处可表示为：

22

1121(,)(,)

(,)(,)

(,)(,)22(,)(0,1,2,,0,1,2,)k j k j k j k j k j k j t x k j u x t u x t u x t u x t a

h

ah

u x t u x h t R x t k j τ

τθτθ++--+''''=+++==±±= (16.2.5) 其中：(,0)()(0,1,2,)k k u x x k ?==±±。由于当,h τ足够小时，在式中略去(,)k j R x t ，就得到一个与方程相近似的差分方程：

1,1,,,0k j k j

k j k j

u u u u a

h

τ

++--+= (16.2.6)

此处，,k j u 可看作是问题的解在节点(,)k j x t 处的近似值。同初值条件：

,0()

(0,1,2,

)k k u x k ?==±± (16.2.7)

结合，就得到求问题的数值解的差分格式。式：

22

21(,)(,)(,)()22

k j k j k j t x ah

R x t u x t u x h t O h τθτθτ''''=+++=+ (16.2.8) 称为差分方程的截断误差。如果一个差分方程的截断误差为()q p R O h τ=+，则称差分方程对

t 是q 阶精度，对x 是p 阶精度的。显然，截断误差的阶数越大，差分方程对微分方程的逼近

越好。

若网格步长趋于0时，差分方程的截断误差也趋于0，则称差分方程与相应的微分方程是相容的。这是用差分方法求解偏微分方程问题的必要条件。

如果当网格步长趋于0时，差分格式的解收敛到相应微分方程定解问题的解，则称这种差分格式是收敛的。

16.3 双曲型方程的差分解法

16.3.1 一阶双曲型方程的差分格式

考虑一阶双曲型方程的初值问题：

00,(,0)()u

u a t x t

x u x x x ????+=>-∞<<+∞?????=-∞<<+∞

?

(16.3.1)

将x t -平面剖分成矩形网格，取x 方向步长为,h t 方向步长为τ，网格线为： (0,1,2,)k x x kh k ===±±, (0,1,2)j t t j j τ=== 为简便，记：

(,)(,)k j k j x t =, (,)(,),()k j k k u k j u x t x ??== 以不同的差商近似偏导数，可以得到方程的不同的差分近似

,1,1,,0k j k j

k j k j

u u u u a

h τ

++--+= (16.3.2)

,1,,1,0k j k j

k j k j

u u u u a h

τ

+---+= (16.3.3)

,1,1,1,02k j k j

k j k j

u u u u a

h

τ

++---+= (16.3.4)

截断误差分别为()O h τ+，()O h τ+与2()O h τ+。结合离散化的初始条件，可以得到几种简单的差分格式：

,1,1,,,0()

(0,1,2,

,0,1,2,)k j k j k j k j k k u u ar u u k j u ?++=--??=±±=?

=?? (16.3.5)

,1,,1,,0()

(0,1,2,,0,1,2,)k j k j k j k j k k

u u ar u u k j u ?+-=--??=±±=?

=?? (16.3.6)

,1,1,1,,0()

(0,1,2,,0,1,2,)2k j k j k j k j k k ar u u u u k j u ?++-?

=--?=±±=?

?=?

(16.3.7)

其中：r h

τ

=

。如果已知第j 层节点上的值,k j u ，按上面三种格式就可求出第1j +层上的值,1k j u +。

因此，这三种格式都是显式格式。

如果对u t ??采用向后差商，u

x

??采用向前差商，则方程可化成：

(,)(,1)

(1,)(,)

()0u k j u k j u k j u k j a

O h h

ττ

--+-+++= (16.3.8)

相应的差分格式为：

,1,1,1,1,0()(0,1,2,

,0,1,2,)k j k j k j k j k k

u u ar u u k j u ?++++=--??=±±=?=?? (16.3.9)

此差分格式是一种隐式格式，必须通过解方程组才能由第j 层节点上的值,k j u 求出第1j +层节点上的值,1k j u +。

例16.3.1 对初值问题：00,(,0)()u u

t x t x

u x x x ????+

=>-∞<<+∞?????=-∞<<+∞

? 其中：101()0200

x x x x ?

==??>??，

用差分格式求其数值解,(1,2,3,4)k j u j =，取1

2

r h

τ

=

=

。 解：记(0,1,2,)k x kh k ==±±，由初始条件：

11,2,1()0

201,2,

k k k x k k ??=--???

===??=?

?

按差分格式：

,1,1,,,0()

(0,1,2,

,0,1,2,)k j k j k j k j k k

u u ar u u k j u ?++=--??=±±=?

=??

计算公式为：,1,1,31

22

k j k j k j u u u ++=-。计算结果略。

如果用差分格式：

,1,,1,,0()

(0,1,2,,0,1,2,)k j k j k j k j k k

u u ar u u k j u ?+-=--??=±±=?

=??

求解，计算公式为：

,1,1,1

()2

k j k j k j u u u +-=+

计算结果略。与准确解：

11(,)200

x t u x t x t x

==??>??

比较知，按前一个差分格式所求得的数值解不收敛到初值问题的解，而后一个差分格式的解收敛到原问题的解。

16.3.2 波动方程的差分格式

对二阶波动方程：

22

222

u u a t x ??=?? (16.3.1) 如果令1u v t ?=

?，2u

v x

?=?，则方程可化成一阶线性双曲型方程组： 21221v v a t x v v t

x ???=????

?

???=???? (16.3.2) 记T 12v (,)v v =，则方程组可表成矩阵形式：

20v v v

10a A t x x ?????==??????? (16.3.3) 矩阵A 有两个不同的特征值a λ=±，故存在非奇异矩阵P ，使得：

100a PAP a -??

==Λ??

-?? (16.3.4) 作变换12w v (,)T P w w ==，方程组可化为：

w w

t x

??

=Λ

??

(16.3.5) 方程组由二个独立的一阶双曲型方程联立而成。因此本节主要讨论一阶双曲型方程的差分解法。

下面给出如下波动方程和边界条件的差分格式：

2

(,)(,)0,0

tt xx

u x t c u x y x a t b

=<<<<(16.3.6) (0,)0,(,)00

(,0)()0

(,0)()0

t

u t u a t t b

u x f t x a

u x g x x a

==≤≤

?

?

=≤≤

?

?=<<

?

(16.3.7)

将矩形{(,):0,0}

R x t x a t b

=≤≤≤≤划分成(1)(1)

n m

-?-个小矩形，长宽分别为：x h

?=，t k

?=，形成一个网格，如图16.3.1所示。

图16.3.1 网格图

可通过差分方程的方法计算近似值：

,

{:1,2,,}

i j

u i n

=在连续行内，2,3,,

j m

=

网格点的真实值为(,)

i i

u x t。

求(,)

tt

u x t和(,)

xx

u x t的中心差分为：

2

2

(,)2(,)(,)

(,)()

tt

u x t k u x t u x t k

u x t O k

k

+-+-

=+(16.3.8)

2

2

(,)2(,)(,)

(,)()

xx

u x h t u x t u x h t

u x t O h

h

+-+-

=+(16.3.9)

在每一行的网格间距是均匀的：

1

i i

x x h

+

=+，而且

1

i i

x x h

-

=-。同时，它在每一列也是均

匀的，1i i t t k +=+，1i i t t k -=-。接下来，将2()O k 和2()O h 去掉，用(16.3.8)和(16.3.9)中的,i j u 近似(,)i j u x t ，并按顺序代入(16.3.6)式，可得到差分方程：

,1,,1

1,,1,2

2

2

22i j i j i j i j i j i j

u u u u u u c k h +-+--+-+= (16.3.10)

可用它来近似(16.3.6)式。为方便起见，可将/r ck h =代入上式，可得：

2,1,,11,,1,2(2)i j i j i j i j i j i j u u u r u u u +-+--+=-+ (16.3.11) 设行j 和1j -的近似值已知，可用上式求网格的行1j +：

22,1,1,1,,1(22)()i j i j i j i j i j u r u r u u u ++--=-++- (16.3.12) 式中，2,3,,1i n =-。根据上式左边的4个已知值可得到近似值,1i j u +，如图16.3.2所示。

图16.3.2 波动方程的空格样板

用上式时，必须注意，如果计算的某个阶段带来的误差最终会越来越小，则方法是稳定的。为了保证上式的稳定性，必然使/1r ck h =≤。还存在其他一些差分方程方法以，称为隐格式法，它们更难实现，但对r 无限制。

16.3.3 差分方法求解波动方程的MATLAB 程序

求解区间{(,):0,0}R x t x a t b =≤≤≤≤，以(16.3.7)为边界条件的波动方程的差分方法程序。

**********************************************************

function U = finedif(f,g,a,b,c,n,m)

%Input - f=u(x,0) as a string 'f'

% - g=ut(x,0) as a string 'g'

% - a and b right endpoints of [0,a] and [0,b]

% - c the constant in the wave equation

% - n and m number of grid points over [0,a] and [0,b]

%Output - U solution matrix; analogous to Table 10.1

%Initialize parameters and U

h = a/(n-1);

k = b/(m-1);

r = c*k/h;

r2=r^2;

r22=r^2/2;

s1 = 1 - r^2;

s2 = 2 - 2*r^2;

U = zeros(n,m);

%Comput first and second rows

for i=2:n-1

U(i,1)=feval(f,h*(i-1));

U(i,2)=s1*feval(f,h*(i-1))+k*feval(g,h*(i-1)) ...

+r22*(feval(f,h*i)+feval(f,h*(i-2)));

end

%Compute remaining rows of U

for j=3:m,

for i=2:(n-1),

U(i,j) = s2*U(i,j-1)+r2*(U(i-1,j-1)+U(i+1,j-1))-U(i,j-2); end

end

U=U';

*******************************************************************

16.4 抛物型方程的差分解法

以一维热传导方程：

220

(0)u u

a a t x

??-=>?? (16.4.1)

为基本模型讨论适用于抛物型方程定解问题的几种差分格式。

16.4.1 差分格式的建立

首先对x t -平面进行网格剖分。分别取,h τ为x 方向与t 方向的步长，用两族平行直线

(0,1,2,)k x x kh k ===±±，(0,1,2)j t t j j τ===，将x t -平面剖分成矩形网格，节点为：(,)(0,1,2,

,0,1,2)k j x t k j =±±=。为简便，记：

(,)(,)k j k j x t =，(,)(,)k j u k j u x t =，()k k x ??=，()k k x ??=，22()j j g g t =，11()j j t λλ=，

22()j j t λλ=。

16.4.2 微分方程的差分近似

在网格内点(,)k j 处，对

u

t

??分别采用向前、向后及中心差商公式：

(,)

(,)

2(,)

(,1)(,)

()(,)(,1)

()

(,1)(,1)

()

2k j k j k j u u k j u k j O t

u u k j u k j O t u u k j u k j O t

ττ

ττ

ττ

?+-=

+??--=

+??+--=

+? (16.4.2)

一维热传导方程可分别表示为：

22

(,1)(,)

(1,)2(,)(1,)

()u k j u k j u k j u k j u k j a

O h h

ττ+-+-+--=+ (16.4.3) 2

2

(,1)(,1)(1,)2(,)(1,)()u k j u k j u k j u k j u k j a O h h

ττ+--+-+--=+ (16.4.4) 22

2

(,1)(,1)(1,)2(,)(1,)()2u k j u k j u k j u k j u k j a O h h

ττ+--+-+--=+ (16.4.5) 由此得到一维热传导方程的不同差分近似：

,1,1,,1,2

20k j k j

k j k j k j

u u u u u a

h τ

++---+-= (16.4.6)

,,1

1,,1,2

20k j k j k j k j k j

u u u u u a

h

τ

-+---+-= (16.4.7)

,1,1

1,,1,2

20k j k j k j k j k j

u u u u u a

h τ

+-+---+-= (16.4.8)

上述差分方程所用到的节点各不相同。其截断误差分别为2()O h τ+，2()O h τ+和

22()O h τ+。因此，它们都与一维热传导方程相容。

如果将式：

22

2

(,1)(,1)(1,)2(,)(1,)()2u k j u k j u k j u k j u k j a O h h

ττ+--+-+--=+中的,k j u 用,1,11

()2

k j k j u u +-+代替，则可得到又一种差分近似：

,1,1

1,,1,11,2

02k j k j k j k j k j k j

u u u u u u a

h

τ

+-++-----+-= (16.4.9)

差分方程用到四个节点。由Taylor 公式容易得出：

2,,1,11

()()2

k j k j k j u u u O τ+-=++ (16.4.10)

故其的截断误差为22

2

2()O h O h ττ??

++ ???

。因而不是对任意的,0h τ→，此差分方程都能逼近热

传导方程：220,(0)u u

a a t x

??-=>??，仅当()o h τ=时，才成立。

综上可知，用不同的差商公式可以得到微分方程的不同的差分近似。构造差分格式的关键在于使其具有相容性、收敛性和稳定性。前面三个方程都具有相容性，而此方程则要在一定条件下才具有相容性。

16.4.3 初、边值条件的处理

为用差分方法求解定解问题初值问题：

220

0,(,0)()u u a t x t x

u x x x ????-=>-∞<<+∞?????=-∞<<+∞

?

(16.4.11)

初边值问题：

2212

00,0(,0)()0(0,)(),(,)()0u u

a t T x l t x u x x x l u t g t u l t g t t T

????-=<<<

=≤≤??==≤≤??? (16.4.12)

还需对定解条件进行离散化。

对初始条件及第一类边界条件，可直接得到： ,0(,0)k k k u u x ?==， (0,1,

0,1,

,)k k n =±=或

0,1,2(0,)(,)j j j n j j j

u u t g u u l t g ====，(0,1,,1)j m =-

式中：,l T

n m h τ

==。

对第二、三类边界条件：

10

122()()

()()

x x l

u t u g t x u t u g t x λλ==???-=?????

??

?+=?????

0t T ≤≤ (16.4.13)

需用差分近似。下面介绍两种较简单的处理方法。

在左边界(0)x =处用向前差商近似偏导数u x ??，在右边界()x l =处用向后差商近似u

x

??，即：

(0,)(,)(1,)(0,)

()

(,)(1,)

()

j n j u u j u j O h x h u u n j u n j O h x h

?-=+??--=+?，(0,1,,)j m = (16.4.14)

则得边界条件的差分近似为：

1,0,10,1,1,2,2j j

j j j n j n j j n j j

u u u g h

u u u g

h λλ--?-=???-?+=??，(0,1,

,)j m = (16.4.15)

其截断误差为()O h 。

(2) 用中心差商近似

u

x

??,即：

2(0,)2(,)(1,)(1,)

()

2(1,)(1,)

()

2j n j u u j u j O h x h u u n j u n j O h x h

?--=+??+--=+?，(0,1,,)j m = (16.4.16)

则得边界条件的差分近似为：

1,1,10,11,1,2,22j j

j j j n j n j j n j j

u u u g h

u u u g

h λλ-+--?-=???-?+=??

， (0,1,

,)j m = (16.4.17)

其截断误差为2()O h 。误差的阶数提高了，但出现定解区域外的节点(1,)j -和(1,)n j +，这就需要将解拓展到定解区域外。可以通过用内节点上的u 值插值求出1,j u -和1,n j u +，也可以假定热传导方程在边界上也成立，将差分方程扩展到边界节点上，由此消去j u ,1-和j n u ,1+。

16.4.4 几种常用的差分格式

以热传导方程的初边值问题：

2212

00,0(,0)()0(0,)(),(,)()0u u

a t T x l t x u x x x l u t g t u l t g t t T

????-=<<<

=≤≤??==≤≤??? (16.4.18)

为例给出几种常用的差分格式。

(1) 古典显式格式 令2

a r h τ

=，则：

,1,1,,1,2

20k j k j

k j k j k j

u u u u u a

h τ

++---+-= (16.4.19)

可改写成：,11,,1,(12)k j k j k j k j u ru r u ru ++-=+-+ 。将其与初始条件及第一类边界条件： ,0(,0)k k k u u x ?==，(0,1,

0,1,

,)k k n =±=或

0,1,2(0,)(,)j j j n j j j

u u t g u u l t g ====，(0,1,,1)j m =-

结合，我们得到求解此问题的一种差分格式：

,11,,1,,00,1,2(12)(1,2,,1,0,1,,1

(0,1,,),(0,1,

,)

k j k j k j k j k k

j

j n j j u ru r u ru k n j m u k n u g u g

j m ?++-?=+-+=-=-?

==??===? (16.4.20)

由于第0层(0)j =上节点处的u 值已知)(0,k k u ?=，由此即可算出u 在第一层(1)j =上节点处的近似值,1k u 。重复使用此式，可以逐层计算出所有的,k j u ，因此此差分格式称为典显式格式。又因式中只出现相邻两个时间层的节点，故此式是二层显式格式。

(2) 古典隐式格式 将式：

,,1

1,,1,2

20k j k j k j k j k j

u u u u u a

h

τ

-+---+-= (16.4.21)

整理并与初始条件及第一类边界条件式联立，得差分格式如下：

,1,1,1,11,1,00,1,2(2)(1,2,

,1,0,1,,1)

(0,1,,)

,(0,1,,)k j k j k j k j k j k k j j n j j u u r u u u k n j m u k n u g u g j m ?++++-+?=+-+=-=-?

==??===? (16.4.22)

其中：2a r h

τ=

。虽然第0层上的u 值仍为已知，但不能由上式直接计算以上各层节点上的值,k j u ，必须通过解下列线性方程组：

1,11,112,12,2,12,1,11,2112012001212j j j j j n j n j n j n j j r r u u rg r r r u u u u u u rg r r r r r +++-+--+-++-??

??+-+-??????????????????????=??????????????????+-+-????????-+??

(16.4.23) 式中，(0,1,

,1)j m =-。 才能由,k j u 计算,1k j u +，故此差分格式称为古典隐式格式。此方程组

是三对角方程组，且系数矩阵严格对角占优，故解存在唯一。

(3) Richardson 格式 Richardson 格式是将式：

,1,1

1,,1,2

20k j k j k j k j k j

u u u u u a

h τ

+-+---+-=

整理后与初始条件及第一类边界条件式联立。其计算公式为：

,1,11,,1,,00,1,22(2)(1,2,,1,1,2,,1

(0,1,,)

,(0,1,,)k j k j k j k j k j k k j

j n j j u u r u u u k n j m u k n u g u g j m ?+-+-?=+-+=-=-?

==??===? 6.4.24)

这种差分格式中所涉及的节点出现在1,,1+-j j j 三层上，故为三层显式格式。Richardson 格式是一种完全不稳定的差分格式，因此它在实际计算中是不能采用的。

(4) 杜福特-弗兰克尔 (DoFort-Frankel) 格式 DoFort-Frankel 格式也是三层显式格式，它是由式：

,1,1

1,,1,11,2

02k j k j k j k j k j k j

u u u u u u a

h

τ

+-++-----+-=

与初始条件及第一类边界条件式结合得到的。具体形式如下：

,11,1,,1,00,1,2212()(1,2,,1,1,2,,1)

1212(0,1,,)

,(0,1,,)k j k j

k j k j k k j j n j j r r u u u u k n j m r r

u k n u g u g j m ?++---?=++=-=-?++?

==??===??

(16.4.25)

用这种格式求解时，除了第0层上的值,0k u 由初值条件得到，必须先用二层格式求出第1层上的值1,k u ，然后再按上式逐层计算,(2,3,

,)k j u j m =。

(5) 六点隐式格式 对二阶中心差商公式：

222

(,)(1,)2(,)(1,)

k j u u k j u k j u k j x h ?+-+-=?

如果用22u x ??在(,1)k j +与(,)k j 处的二阶中心差商的平均值近似22u

x

??在1(,)2k j +处的值，即：

21,1,11,11,,1,222

1(,)

2

2,2()2k j k j k j k j k j k j

k j u u u u u u u

O h x h +++-++-+-++-+?=

+?

同时

u t ??在点)2

1

,(+j k 处的值也用中心差商：

(,)(,1)(,1)

()2k j u u k j u k j O t ττ

?+--=+? 近似，即：

2,1,2

1(,)

2

k j k j

k j u u u

x τ

++-?=

?

这样又得到热传导方程的一种差分近似：

,1,1,1,11,11,,1,2(2,2)02k j k j

k j k j k j k j k j k j u u a

u u u u u u h

τ

++++-++---

-++-+= (16.4.26) 其截断误差为)(22h O +τ，将上式与初始条件及第一类边界条件式联立并整理，得差分格式：

,1

,1,1,11,11,,1,,00,1,2(2)(2)22

(1,2,,1,0,1,,1)(0,1,,),(0,1,,)k j k j k j k j k j k j k j k j k k j j n j j

r r u u u u u u u u k n j m u k n u g u g j m ?++++-++-?

=+-++-+???=-=-??==?

===?? (14.5.27) 此格式涉及到六个节点，它又是隐式格式，故称为六点隐式格式。与古典隐式格式类似，用六点格式由第j 层的值,k j u 计算第1j +层的值,1k j u +时，需求解三对角方程组：

1,12,12,11,11,0,12,1,0,2,3,2,1,,2,1,2,3,1/20/21/200/21/2/21(2)22(2)2

(2)2j j n j n j j j j j j j j j j n j n j n j n j n r r u r r r u u u r r r r r r r u u u u u r u u u u r u u u u u ++-+-++----+-??

??-+-?????????????????????????

?-+-????

??-+??

++-++-+=+-+1,,1,1,2,(2)22j n j n j n j n j r r u u u u +--??????

????????

????

????++-+???? (16.4.28) 式中，(0,1,

,1)j m =-。此方程组的系数矩阵严格对角占优，故仍可用追赶法求解。

例16.4.1 用古典显式格式求初边值问题。

2220

03,03(,0)03

(0,)0,(3,)903u u

t x t x u x x x u t u t t ???-=<<<

=≤≤??==≤≤???

的数值解，取1,0.5h τ==。

解：这里：21,0.5a a r h

τ

===，2()x x ?=，12()0,()9g t g t ==。 由格式：

,11,,1,,00,1,2(12)(1,2,,1,0,1,,1

(0,1,,),(0,1,

,)

k j k j k j k j k k

j j n j j u ru r u ru k n j m u k n u g u g

j m ?++-?=+-+=-=-?

==??===?

可得到：

,11,1,2

2,00,3,0.5()

(1,2,0,1,,5)(0,1,2,3)0,9

(0,1,

,6)

k j k j k j k k

j j u u u k j u x k k u u j ++-=+==??===??===?

将初值,0k u 代入上式，即可算出：

1,12,00,00.5()0.5(40)2u u u =?+=?+= 2,13,01,00.5()0.5(91)5u u u =?+=?+=

将边界条件0,10u =，3,19u =及上述结果代入又可求得：

0,21,22,23,20, 2.5, 5.5,9u u u u ====

如此逐层计算，得全部节点上的数值解为：

0,31,32,33,30,40, 2.75, 5.75,9,0,

u u u u u ====

=

16.4.5 前向差分法求解热传导方程的MATLAB 程序

设(,0)()u x f x =，其中，0x a ≤≤，而且，12(0,),(,)u t c u a t c ==，其中0t b ≤≤，求解区间{(,):0,0}R x t x a t b =≤≤≤≤内2(,)(,)t xx u x t c u x t =的近似解。

function U=forwdif(f,c1,c2,a,b,c,n,m)

%Input - f=u(x,0) as a string 'f' % - c1=u(0,t) and c2=u(a,t)

% - a and b right endpoints of [0,a] and [0,b] % - c the constant in the heat equation

% - n and m number of grid points over [0,a] and [0,b] %Output - U solution matrix; analogous to Table 10.4

%Initialize parameters and U

h=a/(n-1); k=b/(m-1); r=c^2*k/h^2; s=1-2*r;

U=zeros(n,m);

%Boundary conditions

U(1,1:m)=c1; U(n,1:m)=c2;

%Generate first row

U(2:n-1,1)=feval(f,h:h:(n-2)*h)';

%Generate remaining rows of U

for j=2:m

for i=2:n-1

U(i,j)=s*U(i,j-1)+r*(U(i-1,j-1)+U(i+1,j-1)); end end

U=U';

*****************************************************

16.4.6 Crank-Nicholson 求解热传导方程的MATLAB 程序

设(,0)()u x f x =，其中，0x a ≤≤，而且，12(0,),(,)u t c u a t c ==，其中0t b ≤≤，求解区间{(,):0,0}R x t x a t b =≤≤≤≤内2(,)(,)t xx u x t c u x t =的近似解。

Matlab求解微分方程组及偏微分方程组

第四讲 Matlab 求解微分方程（组） 理论介绍：Matlab 求解微分方程（组）命令 求解实例：Matlab 求解微分方程（组）实例 实际应用问题通过数学建模所归纳得到的方程，绝大多数都是微分方程，真正能得到代数方程的机会很少.另一方面，能够求解的微分方程也是十分有限的，特别是高阶方程和偏微分方程（组）.这就要求我们必须研究微分方程（组）的解法：解析解法和数值解法. 一．相关函数、命令及简介 1.在Matlab 中，用大写字母D 表示导数，Dy 表示y 关于自变量的一阶导数，D2y 表示y 关于自变量的二阶导数，依此类推.函数dsolve 用来解决常微分方程（组）的求解问题，调用格式为： X=dsolve(‘eqn1’,’eqn2’,…) 函数dsolve 用来解符号常微分方程、方程组，如果没有初始条件，则求出通解，如果有初始条件，则求出特解. 注意，系统缺省的自变量为t 2.函数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 用初始条件0y 求解. (3)如果要获得微分方程问题在其他指定时间点012,,,,f t t t t 上的解，则令 tspan 012[,,, ]f t t t t =(要求是单调的). (4)因为没有一种算法可以有效的解决所有的ODE 问题，为此，Matlab 提供了多种求解器solver ，对于不同的ODE 问题，采用不同的solver.

用MATLAB解常微分方程

实验四 求微分方程的解 一、问题背景与实验目的 实际应用问题通过数学建模所归纳而得到的方程，绝大多数都是微分方程，真正能得到代数方程的机会很少．另一方面，能够求解的微分方程也是十分有限的，特别是高阶方程和偏微分方程（组）．这就要求我们必须研究微分方程（组）的解法，既要研究微分方程（组）的解析解法（精确解），更要研究微分方程（组）的数值解法（近似解）． 对微分方程（组）的解析解法(精确解)，Matlab 有专门的函数可以用，本实验将作一定的介绍． 本实验将主要研究微分方程(组)的数值解法（近似解），重点介绍 Euler 折线法． 二、相关函数（命令）及简介 1．dsolve ('equ1','equ2',…)：Matlab 求微分方程的解析解．equ1、equ2、…为方程（或条件）．写方程（或条件）时用 Dy 表示y 关于自变量的一阶导数，用用 D2y 表示 y 关于自变量的二阶导数，依此类推． 2．simplify(s )：对表达式 s 使用 maple 的化简规则进行化简． 例如： syms x simplify(sin(x)^2 + cos(x)^2) ans=1 3．[r,how]=simple(s)：由于 Matlab 提供了多种化简规则，simple 命令就是对表达式 s 用各种规则进行化简，然后用 r 返回最简形式，how 返回形成这种形式所用的规则． 例如： syms x [r,how]=simple(cos(x)^2-sin(x)^2) r = cos(2*x) how = combine 4．[T,Y] = solver(odefun,tspan,y 0) 求微分方程的数值解． 说明： (1) 其中的 solver 为命令 ode45、ode23、ode113、ode15s 、ode23s 、ode23t 、ode23tb 之一． (2) odefun 是显式常微分方程：?????==0 0)() ,(y t y y t f dt dy (3) 在积分区间 tspan =],[0f t t 上，从0t 到f t ，用初始条件0y 求解．

用Matlab解微分方程

用Matlab 软件求解微分方程 1．解析解 （1）一阶微分方程 求21y dx dy +=的通解：dsolve('Dy=1+y^2','x') 求y x dx dy -+=21的通解：dsolve('Dy=1+x^2-y','x') 求?????=+=1 )0(12y y dx dy 的特解：dsolve('Dy=1+y^2',’y(0)=1’,'x') （2）高阶微分方程 求解???-='==-+'+''. 2)2(,2)2(,0)(222πππy y y n x y x y x 其中，21=n ，命令为： dsolve('x^2*D2y+x*Dy+(x^2-0.5^2)*y=0','y(pi/2)=2,Dy(pi/2)=-2/pi','x') 求042=-+'-'''x y y y 的通解，命令为： dsolve('D3y-2*Dy+y-4*x=0','x') 输出为： ans=8+4*x+C1*exp(x)+C2*exp(-1/2*(5^(1/2)+1)*x)+C3*exp(1/2*(5^(1/2)-1)*x) （3）一阶微分方程组 求???+-='+='). (3)(4)(),(4)(3)(x g x f x g x g x f x f 的通解：[f,g]=dsolve('Df=3*f+4*g','Dg=-4*f+3*g','x') 输出为： f =exp(3*x)*(cos(4*x)*C1+sin(4*x)*C2) g =-exp(3*x)*(sin(4*x)*C1-cos(4*x)*C2) 若再加上初始条件1)0(,0)0(==g f ，则求特解： [f,g]=dsolve('Df=3*f+4*g','Dg=-4*f+3*g','f(0)=0,g(0)=1','x') 输出为： f =exp(3*x)*sin(4*x) g =exp(3*x)*cos(4*x) 2．数值解 （1）一阶微分方程

Matlab求解微分方程(组)及偏微分方程(组)

第四讲Matlab求解微分方程(组) 理论介绍:Matlab求解微分方程(组)命令 求解实例:Matlab求解微分方程(组)实例 实际应用问题通过数学建模所归纳得到得方程,绝大多数都就是微分方程,真正能得到代数方程得机会很少、另一方面,能够求解得微分方程也就是十分有限得,特别就是高阶方程与偏微分方程(组)、这就要求我们必须研究微分方程(组)得解法:解析解法与数值解法、 一.相关函数、命令及简介 1、在Matlab中,用大写字母D表示导数,Dy表示y关于自变量得一阶导数,D2y 表示y关于自变量得二阶导数,依此类推、函数dsolve用来解决常微分方程(组)得求解问题,调用格式为: X=dsolve(‘eqn1’,’eqn2’,…) 函数dsolve用来解符号常微分方程、方程组,如果没有初始条件,则求出通解,如果有初始条件,则求出特解、 注意,系统缺省得自变量为t 2、函数dsolve求解得就是常微分方程得精确解法,也称为常微分方程得符号解、但就是,有大量得常微分方程虽然从理论上讲,其解就是存在得,但我们却无法求出其解析解,此时,我们需要寻求方程得数值解,在求常微分方程数值解方 面,MATLAB具有丰富得函数,我们将其统称为solver,其一般格式为: [T,Y]=solver(odefun,tspan,y0) 说明:(1)solver为命令ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb、ode15i之一、 (2)odefun就是显示微分方程在积分区间tspan上从到用初始条件求解、 (3)如果要获得微分方程问题在其她指定时间点上得解,则令tspan(要求就是单调得)、 (4)因为没有一种算法可以有效得解决所有得ODE问题,为此,Matlab提供了多种求解器solver,对于不同得ODE问题,采用不同得solver、 表1 Matlab中文本文件读写函数

用Matlab软件求常微分方程的解或通解

《高等数学》实验报告 实验人员:系(班): 学号: 姓名: 实验地点:电教楼五号机房 实验名称:Matlab 高等数学实验 实验时间:2014-6-3 16:30--18:30 实验名称:用Matlab 软件求常微分方程的解(或通解) 实验目的:熟练掌握Matlab 软件求常微分方程的解(或通解) 实验内容:(给出实验程序与运行结果) 1、求微分方程的特解. 1、?? ?? ?===+-10)0(,6)0(034'2 2y y y dx dy dx y d 程序:>> dsolve('D2y-4*Dy+3*y','y(0)=6,Dy(0)=10','x') ans = 4*exp(x)+2*exp(3*x) 吕梁学院《高等数学》实验报告 情况试高中

2、?? ???===++0)0(,2)0(044'2 2y y y dx dy dx y d 程序:>>dsolve('4*D2y+4*Dy+y','y(0)=2,Dy(0)=0','x') ans = 2*exp(-1/2*x)+exp(-1/2*x)*x 3、?? ???===++15)0(',0)0(029422y y y dx dy dx y d 程序:>>dsolve('D2y+4*Dy+29*y=0','y(0)=9,Dy(0)=15','x') ans = 33/5*exp(-2*x)*sin(5*x)+9*exp(-2*x)*cos(5*x) 4、?? ???===+-3)0(',0)0(013422y y y dx dy dx y d 程序:>>dsolve('D2y-4*dy+13*y=0','y(0)=0','Dy(0)=3','x') ans = 3/13*sin(13^(1/2)*x)*13^(1/2)-4/13*cos(13^(1/2)*x)*dy+4/13*dy 5、?? ???-===--5)0(',0)0(04322y y y dx dy dx y d 程序:>>dsolve('D2y-3*Dy-4*y','y(0)=0,Dy(0)=-5','x') ans = exp(-x)-exp(4*x)

MATLAB求解常微分方程数值解

利用MATLAB求解常微分方程数值解

目录 1. 内容简介 (1) 2. Euler Method(欧拉法)求解 (1) 2.1. 显式Euler法和隐式Euler法 (2) 2.2. 梯形公式和改进Euler法 (3) 2.3. Euler法实用性 (4) 3. Runge-Kutta Method(龙格库塔法)求解 (5) 3.1. Runge-Kutta基本原理 (5) 3.2. MATLAB中使用Runge-Kutta法的函数 (7) 4. 使用MATLAB求解常微分方程 (7) 4.1. 使用ode45函数求解非刚性常微分方程 (8) 4.2. 刚性常微分方程 (9) 5. 总结 (9) 参考文献 (11) 附录 (12) 1. 显式Euler法数值求解 (12) 2. 改进Euler法数值求解 (12) 3. 四阶四级Runge-Kutta法数值求解 (13) 4.使用ode45求解 (14)

1.内容简介 把《高等工程数学》看了一遍，增加对数学内容的了解，对其中数值解法比较感兴趣，这大概是因为在其它各方面的学习和研究中经常会遇到数值解法的问题。理解模型然后列出微分方程，却对着方程无从下手，无法得出精确结果实在是让人难受的一件事情。 实际问题中更多遇到的是利用数值法求解偏微分方程问题，但考虑到先从常微分方程下手更为简单有效率，所以本文只研究常微分方程的数值解法。把一个工程实际问题弄出精确结果远比弄清楚各种细枝末节更有意思，因此文章中不追求非常严格地证明，而是偏向如何利用工具实际求解出常微分方程的数值解，力求将课程上所学的知识真正地运用到实际方程的求解中去，在以后遇到微分方程的时候能够熟练运用MATLAB得到能够在工程上运用的结果。 文中求解过程中用到MATLAB进行数值求解，主要目的是弄清楚各个函数本质上是如何对常微分方程进行求解的，对各种方法进行MATLAB编程求解，并将求得的数值解与精确解对比，其中源程序在附录中。最后考察MATLAB中各个函数的适用范围，当遇到实际工程问题时能够正确地得到问题的数值解。 2.Euler Method(欧拉法)求解 Euler法求解常微分方程主要包括3种形式，即显式Euler法、隐式Euler法、梯形公式法，本节内容分别介绍这3种方法的具体内容，并在最后对3种方法精度进行对比，讨论Euler法的实用性。 本节考虑实际初值问题 使用解析法，对方程两边同乘以得到下式

Matlab求解微分方程(组)及偏微分方程(组)

第四讲 Matlab 求解微分方程（组） 理论介绍：Matlab 求解微分方程（组）命令 求解实例：Matlab 求解微分方程（组）实例 实际应用问题通过数学建模所归纳得到的方程，绝大多数都是微分方程，真正能得到代数方程的机会很少.另一方面，能够求解的微分方程也是十分有限的，特别是高阶方程和偏微分方程（组）.这就要求我们必须研究微分方程（组）的解法：解析解法和数值解法. 一．相关函数、命令及简介 1.在Matlab 中，用大写字母D 表示导数，Dy 表示y 关于自变量的一阶导数，D2y 表示y 关于自变量的二阶导数，依此类推.函数dsolve 用来解决常微分方程（组）的求解问题，调用格式为： X=dsolve(‘eqn1’,’eqn2’,…) 函数dsolve 用来解符号常微分方程、方程组，如果没有初始条件，则求出通解，如果有初始条件，则求出特解. 注意，系统缺省的自变量为t 2.函数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 用初始条件0y 求解. (3)如果要获得微分方程问题在其他指定时间点012,,, ,f t t t t 上的解，则令 tspan 012[,,,]f t t t t =(要求是单调的). (4)因为没有一种算法可以有效的解决所有的ODE 问题，为此，Matlab 提供

(完整版)实验七用matlab求解常微分方程

实验七 用matlab 求解常微分方程 一、实验目的： 1、熟悉常微分方程的求解方法，了解状态方程的概念； 2、能熟练使用dsolve 函数求常微分方程（组）的解析解； 3、能熟练应用ode45\ode15s 函数分别求常微分方程的非刚性、刚性的数值解； 4、掌握绘制相图的方法 二、预备知识： 1．微分方程的概念 未知的函数以及它的某些阶的导数连同自变量都由一已知方程联系在一起的方程称为微分方程。如果未知函数是一元函数，称为常微分方程。常微分方程的一般形式为 0),,",',,()(=n y y y y t F Λ 如果未知函数是多元函数，成为偏微分方程。联系一些未知函数的一组微分方程组称为微分方程组。微分方程中出现的未知函数的导数的最高阶解数称为微分方程的阶。若方程中未知函数及其各阶导数都是一次的，称为线性常微分方程，一般表示为 ) ()(')()(1)1(1)(t b y t a y t a y t a y n n n n =++++--Λ 若上式中的系数 n i t a i ,,2,1),(Λ=均与t 无关，称之为常系数。 2．常微分方程的解析解 有些微分方程可直接通过积分求解.例如,一解常系数常微分方程1+=y dt dy 可化为 dt y dy =+1,两边积分可得通解为 1-=t ce y .其中c 为任意常数.有些常微分方程可用一些技巧,如分离变量法,积分因子法,常数变异法,降阶法等可化为可积分的方程而求得解析解. 线性常微分方程的解满足叠加原理,从而他们的求解可归结为求一个特解和相应齐次微分方程的通解.一阶变系数线性微分方程总可用这一思路求得显式解。高阶线性常系数微分方程可用特征根法求得相应齐次微分方程的基本解，再用常数变异法求特解。 一阶常微分方程与高阶微分方程可以互化，已给一个n 阶方程 ),,",',()1()(-=n n y y y t f y Λ 设 ) 1(21,,',-===n n y y y y y y Λ，可将上式化为一阶方程组 ?????????====-),,,,(''''2113221n n n n y y y t f y y y y y y y ΛΛ 反过来，在许多情况下，一阶微分方程组也可化为高阶方程。所以一阶微分方程组与高阶常微分方程的理论与方法在许多方面是相通的，一阶常系数线性微分方程组也可用特征根法求解。 3．微分方程的数值解法 除常系数线性微分方程可用特征根法求解，少数特殊方程可用初等积分法求解外，大部分微分方程无限世界，应用中主要依靠数值解法。考虑一阶常微分方程初值问题

用Matlab件求常微分方程解(或通解)

用Matlab件求常微分方程解(或通解)

————————————————————————————————作者：————————————————————————————————日期：

《高等数学》实验报告 实验人员:系(班): 学号: 姓名: 实验地点:电教楼五号机房 实验名称:Matlab 高等数学实验 实验时间:2014-6-3 16:30--18:30 实验名称:用Matlab 软件求常微分方程的解(或通解) 实验目的:熟练掌握Matlab 软件求常微分方程的解(或通解) 实验内容:(给出实验程序与运行结果) 一、求微分方程的特解. 1、?? ???===+-10)0(,6)0(034'22y y y dx dy dx y d 程序:>> dsolve('D2y-4*Dy+3*y','y(0)=6,Dy(0)=10','x') ans = 4*exp(x)+2*exp(3*x) 吕梁学院《高等数学》实验报告

2、?? ???===++0)0(,2)0(044'22y y y dx dy dx y d 程序:>>dsolve('4*D2y+4*Dy+y','y(0)=2,Dy(0)=0','x') ans = 2*exp(-1/2*x)+exp(-1/2*x)*x 3、?? ???===++15)0(',0)0(029422y y y dx dy dx y d 程序:>>dsolve('D2y+4*Dy+29*y=0','y(0)=9,Dy(0)=15','x') ans = 33/5*exp(-2*x)*sin(5*x)+9*exp(-2*x)*cos(5*x) 4、?? ???===+-3)0(',0)0(013422y y y dx dy dx y d 程序:>>dsolve('D2y-4*dy+13*y=0','y(0)=0','Dy(0)=3','x') ans = 3/13*sin(13^(1/2)*x)*13^(1/2)-4/13*cos(13^(1/2)*x)*dy+4/13*dy 5、?? ???-===--5)0(',0)0(04322y y y dx dy dx y d 程序:>>dsolve('D2y-3*Dy-4*y','y(0)=0,Dy(0)=-5','x') ans = exp(-x)-exp(4*x)

matlab求解常微分方程.docx

用matlab求解常微分方程 在MATLAB中，由函数dsolve()解决常微分方程(组)的求解问题，其具体格式如下：r 二dsolve('eql,eq2,???字condl,cond2,?.?；V) 匕ql,eq2,???*为微分方程或微分方程组，，condl,cond2,.??；是初始条件或边界条件，P是独立变量，默认的独立变量是讥 函数dsolve用来解符号常微分方程、方程组，如果没有初始条件，则求出通解，如果有初始条件，则求出特解。 dy _1 例1：求解常微分方程莎一石的MATLAB程序为： dsolve(* Dy=l/(x+y) 1r!x1), 注意，系统缺省的自变量为t,因此这里要把自变量写明。其中：Y=lambertw(X)表示函数关系Y*exp(Y)二X。 例2：求解常微分方程E'-y— 0的MATLAB程序为： Y2=dsolve(1y*D2y-Dy A2=01, 1x f) Y2=dsolve(!D2y*y-Dy A2=0 J )

我们看到有两个解，其中一个是常数0。 dx 心 ？ —+ 5x + y = e dt 空_兀_3『= g2f 例3：求常微分方程组〔力 ' 通解的MATLAB 程序为： [X,Y]=dsolve(f Dx+5*x+y=exp(t),Dy-x-3*y=exp(2*t) 1, 111 ) [X,Y]=dsolve(1 Dx+2 *x-Dy=l0 * cos(t),Dx+Dy+2 *y=4 *exp(- 2*t) T ,f x(0)=2,y(0)=0f ,f t T ) 以上这些都是常微分方程的精确解法，也称为常微分方程的符号解。但是，我们知 道，有大量的常微分方程虽然从理论上讲，其解是存在的，但我们却无法求出其解析 解，此吋，我们需要寻求方程的数值解，在求常微分方程数值解方面，MATLAB 具有丰 富的函数，我们将其统称为solver,其一般格式为： i°cosr, 7 =2 例4：求常微分方程组 y = 0 z 通解的MATLAB 程序为:

MATLABEuler法解常微分方程

Euler法解常微分方程 Euler法解常微分方程算法： Step 1 分别取积分上限、积分下限、步长 Step 2计算判断是否成立，成立转到Step 3，否则继续进行Step 4 Step 3 计算 Step 4 Euler法解常微分方程算程序： function euler2(fun,y0,A,h) %fun--y' %y0---初值 %A----x取值范围 %a----x左区间端点值 %b----x右区间端点值 %h----给定步长 x=min(A); b=max(A); y=y0; while x

Step 3 （1）做显性Euler预测 （2）将带入 Step 4计算判断是否成立，成立返回Step 3，否则继续进行Step 5 Step 5 改进Euler法解常微分方程算程序： function gaijineuler2(fun,y0,A,h) %fun--y' %y0---初值 %A----x取值范围 %a----x左区间端点值 %b----x右区间端点值 %h----给定步长 a=min(A); b=max(A); x=a:h:b; y(1)=y0; for i=1:length(x)-1 w1=feval(fun,x(i),y(i)); y(i+1)=y(i)+h*w1; w2=feval(fun,x(i+1),y(i+1)); y(i+1)=y(i)+h*(w1+w2)/2; end x=x' y=y' 例：用改进Euler法计算下列初值问题（取步长h=0.25） 输入：fun=inline('-x*y^2') gaijineuler2(fun,2,[0 5],0.25) 得到： x = 0.2500 0.5000 0.7500 1.0000 1.2500 1.5000 1.7500 2.0000 2.2500 2.5000 2.7500

数学应用软件作业6-用Matlab求解微分方程(组)的解析解和数值解

数学应用软件作业6-用Matlab 求解微分方程(组)的解析解和数值解

注：上机作业文件夹以自己的班级姓名学号命名，文件夹包括如下上机报告和Matlab程序。 上机报告模板如下： 佛山科学技术学院 上机报告 课程名称数学应用软件 上机项目用Matlab求解微分方程(组)的解析解和数值解 专业班级姓名学号 一. 上机目的 1.了解求微分方程（组）的解的知识。 2.学习Matlab中求微分方程的各种解的函数，如 dsolve命令、ode45函数等等，其中注意把方程化为新的方程的形式。 3.掌握用matlab编写程序解决求解微分方程的问 题。 二. 上机内容 1、求高阶线性齐次方程：y’’’-y’’-3y’+2y=0。 2、求常微分方程组

2 210cos,2 24,0 t t t dx dy x t x dt dt dx dy y e y dt dt = - = ? +-== ?? ? ?++== ?? 3、求解 分别用函数ode45和ode15s计算求解，分别画出图形，图形分别标注标题。 4、求解微分方程 ,1 )0( ,1 '= + + - =y t y y 先求解析解，在[0,1]上作图； 再用ode45求数值解(作图的图形用“o”表示)，在同一副图中作图进行比较，用不同的颜色表示。 三. 上机方法与步骤 给出相应的问题分析及求解方法，并写出Matlab 程序，并有上机程序显示截图。 题1：直接用命令dsolve求解出微分方程的通解。 Matlab程序：

dsolve('D3y-D2y-3*Dy+2*y','x') 题2：将微分方程组改写为 5cos2exp(2) 5cos2exp(2) (0)2,(0)0 dx t t x y xt dy t t x y dt x y ? =+--- ? ? ? =-+-+- ? ? == ? ? ? ， 再用命令dsolve求解微分方程的通解。 Matlab程序： 建立timu2.m如下： [x,y]=dsolve('Dx=5*cos(t)+2*exp(-2*t)-x-y','Dy=-5*cos(t)+2*exp(-2*t)+x-y ','x(0)=2,y(0)=0','t') x=simple(x) y=simple(y)

matlab简介(解常微分方程绘制函数图像)

MATLAB简介 MATLAB是矩阵实验室（Matrix Laboratory）的简称，是美国MathWorks公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境。 一、基本功能 MATLAB是由美国mathworks公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式，代表了当今国际科学计算软件的先进水平。 MATLAB和Mathematica、Maple并称为三大数学软件。它在数学类科技应用软件中在数值计算方面首屈一指。MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。 MATLAB的基本数据单位是矩阵，它的指令表达式与数学、工程中常用的形式十分相似，故用MATLAB来解算问题要比用C，FORTRAN等语言完成相同的事情简捷得多，并且MATLAB也吸收了像Maple等软件的优点，使MATLAB成为一个强大的数学软件。在新的版本中也加入了对C，FORTRAN，C++，JAVA的支持。可以直接调用,用户也可以将自己编写的实用程序导入到MATLAB函数库中方便自己以后调用，此外许多的MATLAB爱好者都编写了一些经典的程序，用户可以直接进行下载就可以用。 二、特点 1) 高效的数值计算及符号计算功能，能使用户从繁杂的数学运算分析中解脱出来； 2) 具有完备的图形处理功能,实现计算结果和编程的可视化; 3) 友好的用户界面及接近数学表达式的自然化语言，使学者易于学习和掌握； 4) 功能丰富的应用工具箱(如信号处理工具箱、通信工具箱等) ,为用户提供了大量方便实用的处理工具。 三、优势 1.友好的工作平台编程环境 MATLAB由一系列工具组成。这些工具方便用户使用MATLAB的函数和文件，其中许多工具采用的是图形用户界面。包括MATLAB桌面和命令窗口、历史命令窗口、编辑器和调试器、路径搜索和用于用户浏览帮助、工作空间、文件的浏览器。随着MATLAB的商业化以及软件本身的不断升级，MATLAB的用户界面也越来越精致，更加接近Windows的标准界面，人机交互性更强，操作更简单。 2.强大的科学计算机数据处理能力 MATLAB是一个包含大量计算算法的集合。其拥有600多个工程中要用到的数学运算函数，可以方便的实现用户所需的各种计算功能，可以用它来代替底层编程语言，如C和C++ 。在计算要求相同的情况下，使用MATLAB的编程工作量会大大减少。

用matlab求解常微分方程

实验六 用matlab 求解常微分方程 1．微分方程的概念 未知的函数以及它的某些阶的导数连同自变量都由一已知方程联系在一起的方程称为微分方程。如果未知函数是一元函数，称为常微分方程。常微分方程的一般形式为 0),,",',,()(=n y y y y t F 如果未知函数是多元函数，成为偏微分方程。联系一些未知函数的一组微分方程组称为微分方程组。微分方程中出现的未知函数的导数的最高阶解数称为微分方程的阶。若方程中未知函数及其各阶导数都是一次的，称为线性常微分方程，一般表示为 )()(')()(1)1(1)(t b y t a y t a y t a y n n n n =++++-- 若上式中的系数n i t a i ,,2,1),( =均与t 无关，称之为常系数。 2．常微分方程的解析解 有些微分方程可直接通过积分求解.例如,一解常系数常微分方程1+=y dt dy 可化为 dt y dy =+1,两边积分可得通解为 1-=t ce y .其中c 为任意常数.有些常微分方程可用一些技巧,如分离变量法,积分因子法,常数变异法,降阶法等可化为可积分的方程而求得解析解. 线性常微分方程的解满足叠加原理,从而他们的求解可归结为求一个特解和相应齐次微分方程的通解.一阶变系数线性微分方程总可用这一思路求得显式解。高阶线性常系数微分方程可用特征根法求得相应齐次微分方程的基本解，再用常数变异法求特解。 一阶常微分方程与高阶微分方程可以互化，已给一个n 阶方程 ),,",',()1()(-=n n y y y t f y 设)1(21,,',-===n n y y y y y y ，可将上式化为一阶方程组 ?????????====-),,,,(''''2113221n n n n y y y t f y y y y y y y 反过来，在许多情况下，一阶微分方程组也可化为高阶方程。所以一阶微分方程组与高阶常微分方程的理论与方法在许多方面是相通的，一阶常系数线性微分方程组也可用特征根法求解。 3．微分方程的数值解法 除常系数线性微分方程可用特征根法求解，少数特殊方程可用初等积分法求解外，大部分微分方程无限世界，应用中主要依靠数值解法。考虑一阶常微分方程初值问题 ???=<<=000)()),(,()('y t y t t t t y t f t y f

利用matlab编写S函数求解微分方程

利用matlab编写S函数求解微分方程自动化专业综合设计报告 自动化专业综合设计报告

函数求解微S编写设计题目：利用 matlab 分方程 自动化系统仿真实验室所在 实验室： 郭卫平 指导教师： 律迪迪学生姓名 200990519114 班级文自0921 学号 成绩评定： 自动化专业综合设计报告

一、设计目的 了解使用simulink的扩展工具——S-函数，s函数可以利用matlab的丰富资源，而不仅仅局限于simulink提供的模块，而用c或c++等语言写的s函数还可以实现对硬件端口的操作，还可以操作windows API 等的，它的魅力在于完美结合了simulink 框图简洁明快的特点和编程灵活方便的优点，提供了增强和扩展sinulink能力的强大机制，同时也是使用RTW实现实时仿真的关键。 二、设计要求 求解解微分方程 y'=y-2x/y 自动化专业综合设计报告 y(0)=1 要求利用matlab编写S函数求解 三、设计内容（可加附页） 【步骤1】获取状态空间表达式。

在matlab中输入 dsolve(‘Dy=y-2*x/y','y(0)=1'， 'x') 得到 y=(2*x+1).^(1/2); 【步骤2】建立s函数的m文件。 利用21·用S函数模板文件。 以下是修改之后的模板文件sfuntmpl.m 的内容。 function [sys,x0,str,ts] = sfuntmpl(t,x,u,flag) %SFUNTMPL S-function M-file General template define you can With % M-file S-functions, you own ordinary differential system equations (ODEs), discrete % equations, and/or just about any type of algorithm to be used within a %

数学模型之微分方程及其MATLAB求解

数学模型之微分方程及其MATLAB求解 ---卫星轨迹等经典例题求解分析1. 考虑初值问题画图 y'''?3y ''?y 'y = 0 y(0) = 0 y '(0) =1 y ' '(0) = ?1 2、 3、 【实验步骤与程序】 1. M -文件建立m函数文件

function y=f(t,x) y=[x(2);x(3);9*x(3)^2+x(1)*x(2)]; 求解微分方程，命令如下： x0=[0;1;-1]; [t,y]=ode45(@mm,[0,2.5],x0); plot(y(:,1),y(:,2)); figure(2); plot3(y(:,1),y(:,2),y(:,3))

2、M -文件建立m函数文件 function dx=appollo(t,x) mu=1/82.45; mustar=1-mu; r1=sqrt((x(1)+mu)^2+x(3)^2); r2=sqrt((x(1)-mustar)^2+x(3)^2); dx=[x(2) 2*x(4)+x(1)-mustar*(x(1)+mu)/r1^3-mu*(x(1)-mustar)/r2^3 x(4) -2*x(2)+x(3)-mustar*x(3)/r1^3-mu*x(3)/r2^3];

求解微分方程，命令如下： x0=[1.2;0;0;-1.04935751]; options=odeset('reltol',1e-8); [t,y]=ode45(@appollo,[0,20],x0,options); plot(y(:,1),y(:,3)) title('Appollo卫星运动轨迹') xlabel('x') ylabel('y')

Matlab解微分方程(ODE+PDE)

常微分方程： 1 ODE解算器简介(ode**) 2 微分方程转换 3 刚性/非刚性问题(Stiff/Nonstiff) 4 隐式微分方程(IDE) 5 微分代数方程(DAE) 6 延迟微分方程(DDE) 7 边值问题(BVP) 偏微分方程(PDEs)Matlab解法 偏微分方程： 1 一般偏微分方程组(PDEs)的命令行求解 2 特殊偏微分方程(PDEs)的PDEtool求解 3 陆君安《偏微分方程的MATLAB解法 先来认识下常微分方程(ODE)初值问题解算器(solver) [T,Y,TE,YE,IE] = odesolver(odefun,tspan,y0,options) sxint = deval(sol,xint) Matlab中提供了以下解算器： 输入参数： odefun：微分方程的Matlab语言描述函数，必须是函数句柄或者字符串，必须写成Matlab

规范格式(也就是一阶显示微分方程组)，这个具体在后面讲解 tspan=[t0 tf]或者[t0,t1,…tf]：微分变量的范围，两者都是根据t0和tf的值自动选择步长，只是前者返回所有计算点的微分值，而后者只返回指定的点的微分值，一定要注意对于后者tspan必须严格单调，还有就是两者数据存储时使用的内存不同(明显前者多)，其它没有任何本质的区别 y0=[y(0),y’(0),y’’(0)…]：微分方程初值，依次输入所有状态变量的初值，什么是状态变量在后面有介绍 options：微分优化参数，是一个结构体，使用odeset可以设置具体参数，详细内容查看帮助 输出参数： T：时间列向量，也就是ode**计算微分方程的值的点 Y：二维数组，第i列表示第i个状态变量的值，行数与T一致 在求解ODE时，我们还会用到deval()函数，deval的作用就是通过结构体solution计算t 对应x值，和polyval之类的很相似！ 参数格式如下： sol：就是上次调用ode**函数得道的结构体解 xint：需要计算的点，可以是标量或者向量，但是必须在tspan范围内 该函数的好处就是如果我想知道t=t0时的y值，不需要重新使用ode计算，而直接使用上次计算的得道solution就可以 [教程] 微分方程转换为一阶显示微分方程组方法 好，上面我们把Matlab中的常微分方程(ODE)的解算器讲解的差不多了，下面我们就具体开始介绍如何使用上面的知识吧！ 现实总是残酷的，要得到就必须先付出，不可能所有的ODE一拿来就可以直接使用，因此，在使用ODE解算器之前，我们需要做的第一步，也是最重要的一步，借助状态变量将微分

matlabeuler法解常微分方程

Euler 法解常微分方程 Euler 法解常微分方程算法： Step 1 分别取积分上限、积分下限、步长 Step 2计算h n n +=判断b n ≤是否成立，成立转到Step 3，否则继续进行Step 4 Step 3 计算),(1n n n n y x hf y y +=+ Step 4 ),(1n n n n y x hf y y +=+ Euler 法解常微分方程算程序： function euler2(fun,y0,A,h) %fun--y' %y0---初值 %A----x 取值范围 %a----x 左区间端点值 %b----x 右区间端点值 %h----给定步长 x=min(A); b=max(A); y=y0; while x

指导教师： 年 月 日 改进Euelr 法解常微分方程 改进Euler 法解常微分方程算法： Step 1 分别取积分上限、积分下限、步长 Step 2 取一个以h 为步长，a ，b 分别为左右端点的矩阵 Step 3 （1）做显性Euler 预测),(1n n i i y x hf y y +=+ （2）将1+i y 带入)],(),([2 h 111+++++=i i i i i i y x f y x f y y Step 4计算h n n +=判断b n ≤是否成立，成立返回Step 3，否则继续进行Step 5 Step 5 )],(),([2 h 111+++++=i i i i i i y x f y x f y y 改进Euler 法解常微分方程算程序： function gaijineuler2(fun,y0,A,h) %fun--y' %y0---初值 %A----x 取值范围 %a----x 左区间端点值 %b----x 右区间端点值 %h----给定步长 a=min(A); b=max(A); x=a:h:b; y(1)=y0; for i=1:length(x)-1 w1=feval(fun,x(i),y(i)); y(i+1)=y(i)+h*w1; w2=feval(fun,x(i+1),y(i+1)); y(i+1)=y(i)+h*(w1+w2)/2; end x=x'

Matlab求解常微分方程边值问题的方法

Matlab 求解常微分方程边值问题的方法：bvp4c 函数 常微分方程的边值问题，即boundary value problems ，简称BVP 问题，是指表达形式为 (,)((),())0'=??=?y f x y g y a y b 或(,,)((),(),)0'=??=? y f x y p g y a y b p 的方程组（p 是未知参数），在MA TLAB 中使用积分器bvp4c 来求解。 [命令函数] bvp4c [调用格式] sol=bvp4c(odefun,bcfun,solinit,options,p1,p2,…) sol 为一结构体，sol.x 、sol.y 、sol.yp 分别是所选择的网格点及其对应的y(x)与y'(x)数值; bvp4c 为带边值条件常微分方程积分器的函数命令；odefun 为描述微分方程组的函数文件；bcfun 为计算边界条件g(f(a),f(b),p)=0的函数文件；solinit 为一结构体，solinit.x 与solinit.y 分别是初始网格的有序节点与初始估计值，边界值条件分别对应a=solinit.x(l)和b=solinit.x(end)； options 为bvpset 命令设定的可选函数，可采用系统默认值；p1, p2…为未知参数。 例 求常微分方程0''+=y y 在(0)2=y 与(4)2=-y 时的数值解。 [解题过程] 仍使用常用方法改变方程的形式： 令1=y y ，21'=y y ，则原方程等价于标准形式的方程组1221 ?'=??'=-??y y y y ； 将其写为函数文件twoode.m ； 同时写出边界条件函数对应文件twobc.m ； 分别使用结构solinit 和命令bvp4c 确定y-x 的关系； 作出y-x 的关系曲线图。 [算例代码] solinit =bvpinit(linspace(0,4,5),[1 0]); % linspace(0,4,5)为初始网格，[1,0]为初始估计值 sol=bvp4c(@twoode,@twobc,solinit); % twoode 与twobc 分别为微分方程与边界条件的函数，solinit 为结构 x=linspace(0,4); %确定x 范围 y=deval(sol,x); %确定y 范围 plot(x,y(1,:)); %画出y-x 的图形 %定义twoode 函数（下述代码另存为工作目录下的twoode.m 文件） function dydx= twoode(x,y) %微分方程函数的定义 dydx =[y(2) -abs(y(1))]; %定义twobc 函数（下述代码另存为工作目录下的twobc.m 文件） function res= twobc(ya,yb); %边界条件函数的定义 res=[ya(1);yb(1)+2];