常微分方程计算

计 算 题（每题10分）1、求解微分方程2'22x y xy xe -+=。

2、试用逐次逼近法求方程2y x dxdy+=通过点(0,0)的第三次近似解. 3、求解方程'2x y y y e -''+-=的通解4、求方程组dx dt ydydtx y ==+⎧⎨⎪⎩⎪2的通解5、求解微分方程'24y xy x +=6、试用逐次逼近法求方程2y x dxdy-=通过点(1,0)的第二次近似解。

7、求解方程''+-=-y y y e x '22的通解8、求方程组dxdt x ydydtx y =+=+⎧⎨⎪⎩⎪234的通解9、求解微分方程xy y x '-2=24 10、试用逐次逼近法求方程2y x dxdy-=通过(0,0)的第三次近似解. 11、求解方程''+-=-y y y e x '24的通解12、求方程组dxdtx y dydtx y =+=+⎧⎨⎪⎩⎪2332的通解13、求解微分方程x y y e x (')-=14、试用逐次逼近法求方程22x y dxdy+=通过点(0,0)的第三次逼近解. 15、求解方程''+-=--y y y e x '22的通解16、求解方程x e y y y -=-+''32 的通解17、求方程组⎪⎩⎪⎨⎧-+=-+=yx dt dydtdx x y dt dy dt dx243452的通解 18、解微分方程22(1)(1)0x y dx y x dy -+-= 19、试用逐次逼近法求方程2dyx y dx=-满足初始条件(0)0y =的近似解：0123(),(),(),()x x x x ϕϕϕϕ.20、利用逐次逼近法，求方程22dyy x dx=-适合初值条件(0)1y =的近似解：012(),(),()x x x ϕϕϕ。

常微分方程的simpson公式1. 概述常微分方程是描述自变量的一阶导数和未知函数之间关系的方程。

在实际问题中，常微分方程的求解十分重要，可以描述种种自然现象的演化规律，也是物理、生物、经济等多个学科的基础。

在数值求解常微分方程时，simpson公式是一种被广泛应用的数值积分方法，能够在离散的点上对函数进行数值积分。

2. Simpson公式的原理和应用Simpson公式是一种数值积分方法，其原理基于积分中值定理和牛顿-卡索公式。

对于一个给定的函数 f(x)，我们希望计算其在区间 [a, b] 上的定积分∫(a到b) f(x)dx。

Simpson公式给出的积分近似值为\[ \int_{a}^{b} f(x)dx \approx \frac{h}{3} [f(x_0) + 4f(x_1) + 2f(x_2) + \cdots + 2f(x_{n-2}) + 4f(x_{n-1}) + f(x_n)] \]其中，h 是区间 [a, b] 的步长，即 h=(b-a)/n，n 为等分数。

Simpson公式的应用非常广泛，尤其在解常微分方程时扮演重要角色。

利用数值积分方法，我们可以将常微分方程的求解转化为对一阶导数的数值积分。

具体来说，若已知常微分方程的初值条件，我们可以通过适当的步长和积分方法得到未知函数在一系列离散点的近似值，从而得到整个函数的近似解。

3. Simpson公式的优缺点作为一种数值积分方法，Simpson公式具有一定的优缺点。

它的计算简单直观，易于理解和编程实现，因此在实际应用中非常方便。

与其他数值积分方法相比，Simpson公式具有较高的精度和稳定性，特别是对于光滑的函数和较小的步长。

但是，Simpson公式也存在一些不足之处，例如在处理非光滑函数、函数不连续的情况下可能存在较大误差，需要做适当的修正和优化。

4. 结论Simpson公式是一种简单而实用的数值积分方法，广泛应用于常微分方程的求解和其他科学计算领域。

18—1 常微分方程数值解法2§1 引言§2 Euler 方法§3 Runge －Kutta 方法§4 单步法的收敛性与稳定性§5 线性多步法§6 方程组与高阶方程的情况§7 边值问题的数值解法3§1 引言微分方程：关于一个未知函数的方程，方程中含有未知函数的（偏）导数，以及自变量等，其中关于未知函数导数的最高次数称为微分方程的阶数.例如：0)()(')()(''=++−x c y x b y x a x y4实际中，很多问题的数学模型都是微分方程. 常微分方程作为微分方程的基本类型之一，在理论研究与工程实际上应用很广泛. 很多问题的数学模型都可以归结为常微分方程. 很多偏微分方程问题，也可以化为常微分方程问题来近似求解.微分方程的应用情况5对于一个常微分方程：'(,) ,[,]dy y f x y x a b dx==∈为了使解存在，一般要对函数f 施加限制条件，例如要求f 对y 满足Lipschitz 条件：1212(,)(,)f x y f x y L y y −≤−6同时，一个有解的微分方程通常会有无穷多个解例如cos() sin(),dyx y x a a R dx=⇒=+∀∈为了使解唯一，需要加入一个限定条件. 通常会在端点出给出，如下面的初值问题：(,),[,]()dyf x y x a b dx y a y ⎧=∈⎪⎨⎪=⎩7常微分方程的解是一个函数，但是，只有极少数特殊的方程才能求解出来，绝大多数是不可解的.并且计算机没有办法对函数进行运算. 一般考虑其近似解法，一种是近似解析法，如逼近法、级数解法等，另一种是本章介绍的数值解法.8§2 Euler 方法92－1 Euler 公式对常微分方程初值问题：⎩⎨⎧==00')(),(y x y y x f y 数值求解的关键在于消除其中的导数项——称为离散化. 利用差商近似逼近微分是离散化的一个基本途径.10现在假设求解节点为),,1,0(m i ih a x i "=+=，其中ma b h −=为步长，这些节点相应的函数值为)(,),(1m x y x y ". 在点n x 处，已知))(,()('n n n x y x f x y =用n x 的向前差商nn n n x x x y x y −−++11)()(近似代替)('n x y ，如§1，则得到所谓的Euler 公式1(,)n n n n y y hf x y +=+——单步、显式格式11Euler 公式的局部截断误差：假设)(n n x y y =情况下，11)(++−n n y x y 称为局部截断误差.'''2311''23()()()()()2()(,()(()))2n n n n n n n n n y x y x y y x hy x h O h y x h y x f x y x h O h ++−=+++−−=+故有)(2)(''211n n n x y h y x y ≈−++. 122－2 后退的Euler 公式同样对常微分方程初值问题，在1+n x 点，已知))(,()(111'+++=n n n x y x f x y ，如果用向后差商hx y x y n n )()(1−+代替)(1'+n x y ，则得到后退的Euler 公式：111(,)n n n n y y hf x y +++=+——单步、隐式格式13相对于以上可以直接计算1+n y 的Euler 公式（显式），上式是隐式公式. 一般来讲，显式容易计算，而隐式具有更好的稳定性.求解上述公式，通常使用迭代法：对于给定的初值)0(1+n y，计算(1)()111(,)(0,1,)k k n n n n y y f x y k ++++=+="， 如果)(1lim k n k y +∞→收敛，则其极限必满足上述后退Euler 公式.14局部截断误差：假设)(n n x y y =，则),()(111++++=n n n n y x hf x y y .由于)]()[,())(,(),(1111111+++++++−+=n n n y n n n n x y y x f x y x f y x f η且''''2111(,())()()()()n n n n n f x y x y x y x hy x O h +++==++15则有'2''31111(,)[()]()()()()n y n n n n n n y hf x y y x y x hy x h y x O h η++++=−++++将此式减去式2'''31()()()()()2n n n n h y x y x hy x y x O h +=+++ 可得，2''311111()(,)[()]()()2n n y n n n n h y x y hf x y x y y x O h η+++++−=−−+16考虑到21111(,)()1(,)y n y n hf x O h hf x ηη++=++−，则有22''3''11()()()()22n n n n h h y x y y x O h y x ++−=−+≈−172－3 梯形公式由于上述两个公式的局部截断误差绝对值相等，符号相反，故求其算术平均得到梯形公式：111[(,)(,)]2n n n n n n hy y f x y f x y +++=++——单步、隐式格式18梯形法同样是隐式公式，可用下列迭代公式求解：(0)1(1)()111(,)[(,)(,)]2n n n n k k n n n n n n y y hf x y h y y f x y f x y +++++⎧=+⎪⎨=++⎪⎩局部截断误差：类似于后退Euler ，可计算出)(12)('''311n n n x y h y x y −≈−++192－4 改进的Euler 公式上述用迭代法求解梯形公式虽然提高了精度，但计算量也很大. 实际上常采用的方法是，用Euler 公式求得初始值（预测），然后迭代法仅施行一次（校正）——改进的Euler 公式：1111(,)[(,)(,)]2n n n n n n n n n n y y f x y hy y f x y f x y ++++⎧=+⎪⎨=++⎪⎩20估计上式中第二式当1+n y 为准确值时的局部截断误差：''11113(3)()()(()[()()])2()12n n n n n n n hy x y y x y x y x y x hy x ++++−=−++≈−212－5 Euler 两步公式如果用中心差商hx y x y n n 2)()(11−+−代替)('n x y ，则得Euler 两步公式112(,)n n n n y y hf x y +−=+——两步、显式格式22假设1−n y 及n y 均为准确值，利用Taylor 展式容易计算Euler 两步公式的局部截断误差为：11113(3)()()(()2(,()))()3n n n n n n n y x y y x y x hf x y x h y x +++−−=−+≈23此式与梯形公式相结合，得到如下的预测－校正公式：111112(,)[(,)(,)]2n n n n n n n n n n y y hf x y hy y f x y f x y −++++⎧=+⎪⎨=++⎪⎩假设第一式中的1−n y 及n y ，以及第二式中的n y 及1+n y 均是准确值，则有，2441)()(1111−≈−−++++n n n n y x y y x y 从而可得以下的事后估计式，111111114()()51()()5n n n n n n n n y x y y y y x y y y ++++++++⎧−≈−−⎪⎪⎨⎪−≈−⎪⎩25可以期望，以上式估计的误差作为计算结果的补偿，可以提高计算精度.以n p 及n c 分别表示第n 步的预测值和校正值，则有以下的“预测－改进－校正－改进”方案（其中在1+n p 与1+n c 尚未计算出来的前提下，以n n c p −代替11++−n n c p ：26预测：'112n n n hy y p +=−+预测的改进：)(5411n n n n c p p m −−=++计算：),(11'1+++=n n n m x f m校正：)(2'1'1++++=n n n n m y hy c校正的改进：)(511111++++−+=n n n n c p c y计算：),(11'1+++=n n n y x f y27例 用Euler 方法求解初值问题2'[0,0.6](0)1y y xy x y ⎧=−−∈⎨=⎩取0.2h =，要求保留六位小数. 解：Euler 迭代格式为2210.2()0.80.2k k k k k k k k y y y x y y x y +=+−−=−因此2821000(0.2)0.80.20.8y y y x y ≈=−= 22111(0.4)0.80.20.6144y y y x y ≈=−=23222(0.6)0.80.20.461321y y y x y ≈=−=29例 用改进的Euler 方法求解初值问题2'sin 0[0,0.6](0)1y y y x x y ⎧++=∈⎨=⎩取0.2h =，求(0.2),(0.4)y y 的近似值，要求保留六位小数.解：改进的Euler 格式为212211110.2(sin )0.2(sin sin )2k k k k k k k k k k k k k y y y y x y y y y x y y x +++++⎧=+−−⎪⎨=+−−−−⎪⎩30即，222110.820.08sin 0.1(0.80.2sin )sin k k k k k k k k y y y x y y x x ++=−−−则有1(0.2)0.807285y y ≈=，2(0.4)0.636650y y ≈=31§3 Runge －Kutta 方法Def.1如果一种方法的局部截断误差为)(1+p h O ，则称该方法具有p 阶精度. 323－2 Runge —Kutta 方法的基本思想上述的Taylor 级数法虽然可得到较高精度的近似公式，但计算导数比较麻烦. 这里介绍不用计算导数的方法.))(,()()()('1h x y h x f h x y hx y x y n n n n n θθθ++=+=−+——平均斜率.33如果粗略地以),(n n y x f 作为平均斜率，则得Euler 公式；如果以221K K +作为平均斜率，其中),(1n n y x f K =，),(112hK y x f K n n +=+，则得改进的Euler 公式.343－3 二阶的Runge －Kutta 方法对点n x 和)10(≤<+=+p ph x x n p n ，用这两点斜率的线性组合近似代替平均斜率，则得计算公式：11122121()(,)(,)n n n n n p n y y h K K K f x y K f x y phK λλ++⎧=++⎪=⎨⎪=+⎩35现确定系数p ,,21λλ，使得公式具有二阶精度. 因为，取n y 为()n y x ，则'1(,)(,())'()n n n n n nK f x y f x y x y x y === 再把2K 在),(n n y x 处展开，有36'21(,)(,)n p n n n n K f x y phK f x ph y phy +=+=++代入可得，'2''31122()()n n n n y y hy ph y O h λλλ+=++++'2(,)(,)(,)()n n x n n y n n n f x y f x y ph f x y phy O h =+⋅+⋅+'2(')(,)()n x y n n y ph f f y x y O h =+⋅+⋅+'''2()n n y ph y O h =+⋅+37相比较二阶Taylor 展开''2'12n n n n y h hy y y ++=+，有，⎪⎩⎪⎨⎧==+211221p λλλ满足此条件的公式称为二阶Runge －Kutta 公式.38可以验证改进的Euler 公式属于二阶Runge －Kutta 公式. 下列变形的Euler 公式也是二阶Runge －Kutta 公式：12121(,)(,)22n n n n n n y y hK K f x y h h K f x y K +⎧⎪=+⎪=⎨⎪⎪=++⎩393－4 三阶Runge －Kutta 公式同二阶Runge －Kutta 公式，考虑三点,,(01)n n p n q x x x p q ++≤≤≤试图用它们的斜率321,,K K K 的线性组合近似代替平均斜率，即有如下形式的公式：1112233121312()(,)(,)(,())n n n n n n n n y y h K K K K f x y K f x ph y phK K f x qh y qh rK sK λλλ+=+++⎧⎪=⎪⎨=++⎪⎪=+++⎩40把32,K K 在),(n n y x 处展开，通过与)(1+n x y 在n x 的直接Taylor 展式比较，可确定系数s r q p ,,,,,,321λλλ，满足下式，从而使得上述公式具有三阶精度，41特别地，2,1,1,21,32,61231=−======s r q p λλλ是其一特例.123232223311213161p q p q pqs r s λλλλλλλλ++=⎧⎪⎪+=⎪⎪⎪+=⎨⎪⎪=⎪⎪+=⎪⎩423－5 四阶Runge －Kutta 公式相同的方法，可以导出下列经典的四阶Runge －Kutta 公式：112341213243(22)6(,)(,)22(,)22(,)n n n n n n n n n n h y y K K K K K f x y h h K f x y K h h K f x y K K f x h y hK +⎧=++++⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪=++⎪⎩43例 用经典四阶Runge —Kutta 方法求解初值问题'83[0,0.4](0)1y y x y =−⎧∈⎨=⎩，取0.2h =，求(0.4)y 的近似值，要求保留六位小数.解：四阶Runge —Kutta 格式为44112341211123122241330.2(22)6(,)830.2(,)83(0.1) 5.6 2.120.2(,)83(0.1) 6.32 2.372(,0.2)83(0.2) 4.208 1.578k k k k k k k k k k k kk k k k ky y K K K K K f x y y K f x y K y K yK f x y K y K y K f x y K y K y ++++⎧=++++⎪⎪==−⎪⎪⎪=+=−+=−⎨⎪⎪=+=−+=−⎪⎪⎪=+=−+=−⎩则10.5494 1.2016k k y y +=+，45故12(0.2) 2.3004,(0.4) 2.4654y y y y ≈=≈=.注：由准确解382()33xy x e −=−可得(0.2) 2.300792,(0.4) 2.465871y y ==46§5 线性多步法基本思想：在计算1+i y 之前，已计算出一系列的近似值i y y ,,1"，如果充分利用这些已知信息，可以期望会获得更高精度的)(1+i x y 的近似值1+i y .基本方法：基于数值积分与基于Taylor 展开的构造方法.475－1 基于数值积分的构造方法对方程),('y x f y =两边从i x 到1+i x 积分，则得∫++=+1),()()(1i ix x i i dxy x f x y x y 设)(x P r 是f (x , y )的插值多项式，由此可得以下的一般形式的计算公式：∫++=+1)(1i ix x r i i dxx P y y 48例 取线性插值))(,())(,()(11111+++++−−+−−=i i i i ii i i i i r x y x f x x x x x y x f x x x x x P ，则得到梯形法：)],(),([2111+++++=i i i i i i y x f y x f hy y495－2 Adams 显式公式在区间],[1+i i x x 上利用r ＋1个数据点),(,),,(),,(11r i r i i i i i f x f x f x −−−−"构造插值多项式)(x P r ，由牛顿后插公式（注意到：j i j i j f f −Δ=∇）j i jrj j i r f j t th x P −=Δ⎟⎟⎠⎞⎜⎜⎝⎛−−=+∑0)1()(其中!)1()1(j j s s s j s +−−=⎟⎟⎠⎞⎜⎜⎝⎛". 50可得10rj i i rj i jj y y h f αΔ+−==+∑——Adams 显式公式其中1(1)j j t dt j α−⎛⎞=−⎜⎟⎝⎠∫，它可写成：∑=−++=rj ji rj i i f h y y 01β515－3 Adams 隐式公式在区间],[1+i i x x 上利用r ＋1个数据点),(,),,(),,(1111+−+−++r i r i i i i i f x f x f x "构造插值多项式)(x P r ，由牛顿后插公式101)1()(+−=+Δ⎟⎟⎠⎞⎜⎜⎝⎛−−=+∑j i jrj ji r f j t th x P 可得*11rj i i rj i j j y y h f α+−+==+Δ∑——Adams 隐式公式52其中01(1)jj t dt j −−⎛⎞α=−⎜⎟⎝⎠∫，它又可写成： *11ri i rj i j j y y h f β+−+==+∑535－4 Adams 预测－校正公式以r ＝3时的Adams 显式与隐式公式为例. 此时，显式公式为)9375955(243211−−−+−+−+=i i i i i i f f f f hy y 利用Taylor 展式，容易计算局部截断误差为)(720251)5(5i x y h . 54)5199(242111−−+++−++=i i i i i i f f f f hy y 同样利用Taylor 展开可得，其局部截断误差为5(5)19()720i h y x −. 隐式公式为55⎪⎩⎪⎨⎧+−++=−+−+=−−+++−−−+)519),(9(24)9375955(24211113211i i i i i i i i i i i i i f f f y x f hy y f f f f h y y 注 利用2－5节的相同作法同样可以构造更精确的计算过程.可构造利用显式预测，隐式校正的计算公式：56§6 方程组与高阶方程的情形6－1 一阶方程组常微分方程初值问题为⎩⎨⎧==00)(),('y x y y x f y 此时T m y y y ),,(1"=，Tm f f f ),,(1"=. 此时上述的一切方法均可使用，只是注意y 与f 此时为向量.576－2 化高阶方程为一阶方程组解下列的m 阶方程()(1)'(1)(1)000000(,,',,)(),'(),,()m m m m y f x y y y y x y y x y yx y −−−⎧=⎨===⎩""令)1(21,,',−===m m y y y y y y "，则有58'12'23'1'12(,,,,)m m m m y y y y y yy f x y y y −⎧=⎪=⎪⎪⎨⎪=⎪⎪=⎩#"初始条件为：)1(00'002001)(,,)(,)(−===m m y x y y x y y x y "。

求常微分方程的数值解一、背景介绍常微分方程（Ordinary Differential Equation，ODE）是描述自然界中变化的数学模型。

常微分方程的解析解往往难以求得，因此需要寻找数值解来近似地描述其行为。

求解常微分方程的数值方法主要有欧拉法、改进欧拉法、龙格-库塔法等。

二、数值方法1. 欧拉法欧拉法是最简单的求解常微分方程的数值方法之一。

它基于导数的定义，将微分方程转化为差分方程，通过迭代计算得到近似解。

欧拉法的公式如下：$$y_{n+1}=y_n+f(t_n,y_n)\Delta t$$其中，$y_n$表示第$n$个时间步长处的函数值，$f(t_n,y_n)$表示在$(t_n,y_n)$处的导数，$\Delta t$表示时间步长。

欧拉法具有易于实现和理解的优点，但精度较低。

2. 改进欧拉法（Heun方法）改进欧拉法又称Heun方法或两步龙格-库塔方法，是对欧拉法进行了精度上提升后得到的一种方法。

它利用两个斜率来近似函数值，并通过加权平均来计算下一个时间步长处的函数值。

改进欧拉法的公式如下：$$k_1=f(t_n,y_n)$$$$k_2=f(t_n+\Delta t,y_n+k_1\Delta t)$$$$y_{n+1}=y_n+\frac{1}{2}(k_1+k_2)\Delta t$$改进欧拉法比欧拉法精度更高，但计算量也更大。

3. 龙格-库塔法（RK4方法）龙格-库塔法是求解常微分方程中最常用的数值方法之一。

它通过计算多个斜率来近似函数值，并通过加权平均来计算下一个时间步长处的函数值。

RK4方法是龙格-库塔法中最常用的一种方法，其公式如下：$$k_1=f(t_n,y_n)$$$$k_2=f(t_n+\frac{\Delta t}{2},y_n+\frac{k_1\Delta t}{2})$$ $$k_3=f(t_n+\frac{\Delta t}{2},y_n+\frac{k_2\Delta t}{2})$$ $$k_4=f(t_n+\Delta t,y_n+k_3\Delta t)$$$$y_{n+1}=y_n+\frac{1}{6}(k_1+2k_2+2k_3+k_4)\Delta t$$三、数值求解步骤对于给定的常微分方程，可以通过以下步骤求解其数值解：1. 确定初值条件：确定$t=0$时刻的函数值$y(0)$。

常微分方程的数值解算法常微分方程的数值解算法是一种对常微分方程进行数值计算的方法，这可以帮助我们更好地理解和研究自然现象和工程问题。

在本文中，我们将介绍一些常用的数值解算法，探讨它们的优缺点和适用范围。

常微分方程（ODE）是描述自然现象和工程问题的重要数学工具。

然而，对于许多ODE解析解是无法求出的，因此我们需要通过数值方法对其进行求解。

常微分方程可以写作：y' = f(t, y)其中，y是函数，f是给定的函数，表示y随t的变化率。

这个方程可以写成初始值问题（IVP）的形式：y'(t) = f(t,y(t))，y(t0) = y0其中，y（t0）=y0是方程的初始条件。

解决IVP问题的典型方法是数值方法。

欧拉方法欧拉方法是最简单的一阶数值方法。

在欧拉方法中，我们从初始条件开始，并在t = t0到t = tn的时间内，用以下公式逐步递推求解：y n+1 = y n + hf (t n, y n)其中，f（t n，y n）是点（t n，y n）处的导数， h = tn - tn-1是时间间隔。

欧拉方法的优点是简单易懂，容易实现。

然而，它的缺点是在整个时间段上的精度不一致。

程度取决于使用的时间间隔。

改进的欧拉方法如果我们使用欧拉方法中每个时间段的中间点而不是起始点来估计下一个时间点，精度就会有所提高。

这个方法叫做改进的欧拉方法（或Heun方法）。

公式为：y n+1 = y n + h½[f(t n, y n)+f(tn+1, yn + h f (tn, yn))]这是一个二阶方法，精度比欧拉方法高，但计算量也大一些。

对于易受噪声干扰的问题，改进的欧拉方法是个很好的选择。

Runge-Kutta方法Runge-Kutta方法是ODE计算的最常用的二阶和高阶数值方法之一。

这个方法对定义域内的每个点都计算一个导数。

显式四阶Runge-Kutta方法（RK4）是最常用的Runge-Kutta方法之一，并已得到大量实践的验证。