随机微分方程数值解法.

微分数值解法分别用欧拉方法、改进的欧拉方法、二阶龙格-库塔方法、三阶龙格-库塔方法、四阶龙格-库塔方法，计算下列初值问题，并与精确解对比。

微分方程的初值问题:⎪⎩⎪⎨⎧==00)(),(y x y y x f dx dy1.欧拉方法：用差分代替微分，化为：⎩⎨⎧+=+=++nhx x y x hf y y n n n n n 011),( 2.改进的欧拉方法：⎪⎩⎪⎨⎧+=++=-++-++),()],(),([21111n n n n n n n n n n y x hf y y y x f y x f h y y 3二阶龙格-库塔方法：⎪⎪⎩⎪⎪⎨⎧++==+=+)2,2(),(k 12121κκκh y h x f y x f h y y n n n n n n 4.三阶龙格-库塔方法：⎪⎪⎪⎩⎪⎪⎪⎨⎧+-++=++==+++=+))2(,()2,2(),()4(612131213211κκκκκκκκκh y h x f h y h x f y x f y y n n n n n n n n 5，四阶龙格-库塔方法：⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎨⎧++=++=++==++++=+),()2,2()2,2(),()22(6342312143211κκκκκκκκκκκh y h x f h y h x f h y h x f y x f h y y n n n n nn n n n n例一：⎩⎨⎧=-+=0)0(22'y x x y y 易得方程精确解为：y=2x根据欧拉方法、改进的欧拉方法、二阶龙格-库塔方法、三阶龙格-库塔方法、四阶龙格-库塔方法，利用c++编程（程序见附录），计算结果如下表：例二:⎩⎨⎧=+-=0)0(x2x 24'y y y 易得精确解为y=2x根据欧拉方法、改进的欧拉方法、二阶龙格-库塔方法、三阶龙格-库塔方法、四阶龙格-库塔方法，利用c++（程序见附录），计算结果如下表：例三:⎩⎨⎧=--=1)0(sin x cos 22'y xy y 易得精确解为y=cosx根据欧拉方法、改进的欧拉方法、二阶龙格-库塔方法、三阶龙格-库塔方法、四阶龙格-库塔方法，利用c++（程序与例一类似，故省略见附录），计算结果如下表：附录：例一程序（1）欧拉方法的c++程序：（2）改进的欧拉方法的c++程序：（3）二阶龙格-库塔方法：（4）三阶龙格-库塔方法：（5）四阶龙格-库塔方法：例二程序：（1）欧拉方法的c++程序：（2）改进的欧拉方法的c++程序：（3）二阶龙格-库塔方法：（4）三阶龙格-库塔方法：（5）四阶龙格-库塔方法例三程序：（1）欧拉方法的c++程序：（2）改进的欧拉方法的c++程序：（3）二阶龙格-库塔方法：（4）三阶龙格-库塔方法：（5）四阶龙格-库塔方法：。

微分方程数值解法微分方程数值解法微分方程数值解法【1】摘要：本文结合数例详细阐述了最基本的解决常微分方程初值问题的数值法，即Euler方法、改进Euler法，并进行了对比，总结了它们各自的优点和缺点，为我们深入探究微分方程的其他解法打下了坚实的基础。

关键词:常微分方程数值解法 Euler方法改进Euler法1、Euler方法由微分方程的相关概念可知，初值问题的解就是一条过点的积分曲线，并且在该曲线上任一点处的切线斜率等于函数的值。

根据数值解法的基本思想，我们取等距节点，其中h为步长，在点处，以为斜率作直线交直线于点。

如果步长比较小，那么所作直线与曲线的偏差不会太大，所以可用的近似值，即：，再从点出发，以为斜率作直线，作为的近似值，即：重复上述步骤，就能逐步求出准确解在各节点处的近似值。

一般地，若为的近似值，则过点以为斜率的直线为：从而的近似值为：此公式就是Euler公式。

因为Euler方法的思想是用折线近似代替曲线，所以Euler方法又称Euler折线法。

Euler方法是初值问题数值解中最简单的一种方法，由于它的精度不高，当步数增多时，由于误差的积累，用Euler方法作出的折线可能会越来越偏离曲线。

举例说明：解： ,精确解为：1.2 -0.96 -1 0.041.4 -0.84 -0.933 0.9331.6 -0.64 -0.8 0.161.8 -0.36 -0.6 0.242.0 0 -0.333 0.332.2 0.44 0 0.44通过上表可以比较明显地看出误差随着计算在积累。

2、改进Euler法方法构造在常微分方程初值问题 ,对其从到进行定积分得:用梯形公式将右端的定积分进行近似计算得：用和来分别代替和得计算格式:这就是改进的Euler法。

解：解得：由于 ,是线形函数可以从隐式格式中解出问题的精确解是误差0.2 2.421403 2.422222 0.000813 0.021400.4 2.891825 2.893827 0.00200 0.051830.6 3.422119 3.425789 0.00367 0.094112.0 10.38906 10.43878 0.04872 1.1973通过比较上表的第四列与第五列就能非常明显看出改进Euler方法精度比Euler方法精度高。

微分⽅程数值解法答案包括基本概念，差分格式的构造、截断误差和稳定性，这些内容是贯穿整个教材的主线。

解答问题关键在过程，能够显⽰出你已经掌握了书上的内容，知道了解题⽅法。

这次考试题⽬的类型：20分的选择题，主要是基本概念的理解，后⾯有五个⼤题，包括差分格式的构造、截断误差和稳定性。

习题⼀1．略2． y y x f -=),(，梯形公式：n n n n n n y hh y y y h y y )121(),(2111+-+=+-=+++，所以0122)1(01])121[()121()121(y hh y h h y h h y hhn h h n n n +--+--+-+=+-+==+-+= ，当0→h 时，x n e y -→。

同理可以证明预报-校正法收敛到微分⽅程的解.3．局部截断误差的推导同欧拉公式;整体截断误差:++++++-++≤1),())(,(11111n nx x n n n n n n n dx y x f x y x f R εε11)(++-++≤n n n y x y Lh R ε，这⾥R R n ≤ ⽽111)(+++-=n n n y x y ε，所以 R Lh n n +=-+εε1)1(，不妨设1()]11111[1111101---++-+-+-≤≤-+-=n n n n Lh Lh Lh R Lh Lh R Lh εεε ]1[2)(02)(00-+≤--x X L x X L eLh R eε4．中点公式的局部截断误差： dx x y x f hx y h x f x y x f yx y n n x x n n n n n n))](,(2)(,2())(,([)(11*1?+++-=-++dx x y x f hx y h x f h x y h x f h x y x y dxx y x f hx y h x f hx y h x f h x y h x f x y x f n n n n x x n n n n n n n x x n n n n n n n n))](,(2)(,2())2(,2([)]2()([))](,(2)(,2())2(,2())2(,2())(,([11++-++++'-'=++-+++++-=??++所以上式为+--+''=?++dx hx x x y e n nx x n n n )2()(11θdx x y x f h x y h x f h x y h x f n n n n x x n n n n))](,(2)(,2())2(,2([1++-++?+ 3218)(LMh h x y Lh e n n ≤+''≤+?中点公式的整体截断误差：dx y x f hy h x f x y x f y x y y x y n n x x n n n n n n n n)],(2,2())(,([)()(111?+++-+-=-++dxy x f hy h x f x y x f h x y h x f x y x f hx y h x f x y x f y x y n n n n n n n n x x n n n n n n n n))],(2,2()))(,(2)(,2()))(,(2)(,2())(,([)(1++-+++++-+-=?+因⽽n n n L h Lh R εεε)21(1+++≤+，R L h Lh n n +++≤-122)21(εε≤≤])21()21(1[2)21(1222222022-+++++++--+++n nL h Lh L h Lh Lh Lh RL h Lh ε )1(00-+≤--x X L x X L e LhR eε 5．略 6．略 7．略8．（1）欧拉法：2.0≤h ；四阶Runge-Kutta ⽅法：278.0≤h （2）欧拉法：3 54≤h ；四阶Runge-Kutta ⽅法：3556.5≤h（3）欧拉法：1≤h ；四阶Runge-Kutta ⽅法：278.0≤h 9．略 10．略习题21．略 2．略 3．略4．差分格式写成矩阵形式为：n n M n M n n n M n M n n e u u u u r t r r r t r r r t r r r t u u u u +?--------= --+-+-++12211221121212121 αβαααβαααβαααβ矩阵的特征值为：)cos(221Mj r r t j πααβλ+-?-=，要使格式稳定，则特征值须满⾜t c j ?+≤1λ，即21≤r α5．利⽤泰勒展式可以得到古典隐式差分格式的截断误差为)(2h t O +?。

微分方程数值解及其应用绪论自然界中的许多事物的运动和变化规律都可以用微分方程来描述，因此对工程和科学技术中的实际问题的研究中,常常需要求解微分方程．但往往只有少数较简单和典型的微分方程可求出其解析解，在大多数情况下, 只能用近似法求解，数值解法是一类重要的近似方法．本文主要讨论一阶常微分方程的初值问题的数值解法，探讨这些算法在处理来自生活实际问题中的应用，并结合MATLAB软件，动手编程予以解决．１微分方程的初值问题[1]1.1 预备知识在对生活实际问题的研究中，通常需要考虑一阶微分方程的初值问题dyf ( x, y)dxy( x0 )y0〔1〕这里 f x, y 是矩形区域 R ： x x0a, y y0 b 上的连续函数．对初值问题〔 1〕需要考虑以下问题：方程是否一定有解呢?假设有解，有多少个解呢?下面给出相关的概念与定理．定义 1Lipschitz 条件[1][ 2]：矩形区域 R ： x x0a, y y0 b 上的连续函数f x, y 假设满足：存在常数 L0 ，使得不等式 f x, y1 f x, y2L y1 y 2对所有x, y1 , x, y2R 都成立，那么称f x, y在 R 上关于 y 满足Lipschitz条件 .定理 1解的存在唯一性定理[1] [3]：设 f 在区域 D x, y a x b, y R 上连续，关于 y 满足Lipschitz条件，那么对任意的 x0a, b , y0R ，常微分方程初值问题〔1〕当 x a, b 时存在唯一的连续解 y x ．该定理保证假设一个函数f x, y 关于 y 满足Lipschitz条件，它所对应的微分方程的初值问题就有唯一解 . 在解的存在唯一性得到保证的前提下，自然要考虑方程的求解问题．求解微分方程虽然有多种解析方法， 但根据工程和科学实践问题所得到的微分方程往往很复杂，在很多情况下不能或很难给出解析解，有时即使能求出形式解，也往往因形式过于复杂或计算量太大而不实用，因此从实际问题中归结出来的微分方程主要依靠数值解法．定义2微分方程数值解：对初值问题〔1〕寻求数值解就是寻求解y x 在一系列离散节点上的近似解y 0 , y 1, y 2 ,, y n , y n 1,，相邻两个节点的间距h n x n 1x n称为步长 . 在一般情况下假定h ih i0,1,为常数，这时节点为x nx 0nh, n0,1,2,．要求微分方程数值解，首先要建立数值算法，即对初值问题〔 1〕中的方程离散化，建立求解数值解法的递推公式．一类是计算 y n 1 时只用到前一点的值 y n ，称为单步法；另一类是用到 y n 1 前面 k 点的值 y n , y n 1, , y n k 1 称为 k 步法 .对初值问题〔 1〕式的单步法可用一般形式表示为yn 1y n h (x n , y n , y n 1 , h) ，其中多元函数 与 f x, y 有关，当含有 y n 1 时，方法是隐式的；假设 中不含 y n 1 ，那么为显式方法，所以显式单步法可表示为yn 1y n h (x n , y n ,h).〔2〕设 y x 是初值问题〔 1〕的准确解，称 T n 1 y x n 1 y x n h( x n , y x n , h) 为显式单步法〔 2〕的局部截断误差 . 假设存在最大正整数 p , 使显式单步法〔 2〕式的局部截断误差满足 T n 1 y x h y xhx, y, h O h p 1 ，那么称〔 2〕式有 p 阶精度 .1.2 几种常用的数值解法及其分析、比拟1.2.1 欧拉法与后退欧拉法1〕欧拉法：欧拉曾简单地用差分代替微分，即利用公式将初值问题〔 1〕离散化，那么问题〔 1〕可化为y n 1 y n h f (x n , y n ), x n x 0 n h ，〔3〕此方法称为欧拉法 .欧拉方法的几何意义在数值计算公式中表达了出来. 在 xy 平面上，一阶微分方程的解y y x 称作它的积分曲线 . 积分曲线上一点x, y 的切线斜率等于函数f x, y ，按函数 f x, y 在 xy 平面上建立一个方向场，那么，积分曲线上每一点的切线方向均与方向场在该点的方向相一致.基于上述几何解释，从初始点 P0 ( x0 , y0 ) 出发，先依方向场在该点的方向上推进到x x1上一点 P1，再从 P1依方向场的方向推进到 x x2上一点 P2，循环前进便作出一条折线 P PP，因此欧拉方法又称为折线法. 假设初值y，那么由〔 3〕式可逐步算0 1 20出为了分析计算公式的精确度，通常可用泰勒展开将y x1在 x n处展开，那么有ny x n 1 y x n h y x nh2y''n,n x n , x n 1 . y x n h2在 y n y x n的前提下， f x n , y n f x n , y x n y x n . 可得欧拉法〔 3〕的误差为容易看出，欧拉法〔 3〕式具有一阶精度 .2〕向后欧拉方法：如果对微分方程〔1〕从x n到x n 1积分，得y x n 1xn 1〔4〕y x n f t , y t dx ，x n如果〔 4〕式右端积分用右矩形公式h f x n 1 , y x n 1近似，那么得到另一个公式yn 1y n hf x n 1 , y n 1，〔5〕称为后退欧拉法 .值得一提的是 : 后退欧拉法与欧拉公式有着本质的区别，后者是关于y n 1的直接计算公式，它是显式的，而〔5〕式的右端含有关于y n 1的表达式，它是隐式的.在利用后退欧拉法时，我们通常利用迭代法求解，实质就是逐步显示化. 具体迭代过程如下：首先利用欧拉公式 y n (0)1 y n h f ( x n , y n ) 给出迭代初值 y n (0)1 ，把它代入〔 5〕式的y(1)y n h f ( x n 1, y (0)右端，使之转化为显式，直接计算得n 1n 1 ) . 如此反复进行，得y n (k 1 1) y n h f ( x n 1 , y n (k 1)) k0,1, ，那么得到后退欧拉法的迭代公式y n (0)1 y n h f ( x n , y n ) ，y n (k 11)y n h f ( x n 1, y n (k 1) ) 可以看出，后退欧拉法具有一阶精度，且计算比拟麻烦 .1.2.2 梯形方法为得到比欧拉法精确度高的计算公式，在等式〔 4〕式右端积分中假设用梯形求积公式近似，并用 y n 代替 y x n ， y n 1 代替 y x n 1 ，那么得y n 1ynh f x n 1, y n 1，f x n , y n2〔6〕称其为梯形方法 .梯形方法与后退欧拉法一样， 都是隐式单步法， 可用迭代法求解， 其迭代公式为y n (0)1 y n h f (x n , y n ) .〔7〕y n (k 11)y n hf x n , y nf x n 1, y n (k 1)2为了分析梯形公式的收敛性，将〔6〕与〔 7〕式相减，得(k 1)h x n 1, yn 1f xn(k ) ， k 0,1,2,yn 1yn 1f1 , y n 12因为 fx, y 满足 Lipschitz 条件，于是有 y n 1 y n(k11)hLyn 1y n (k 1) ，其中 L 为 f x, y2关于 y 的 Lipschitz 常 数 . 如 果选 取 h 充分 小 ， 使 得hL1 ， 那么当k时有2y n (k 11 )y n 1 ，这说明迭代过程〔 7〕式是收敛的 [4] . 容易推导得出梯形法〔7〕式是二阶方法 .经分析，梯形方法虽然提高了精度，但是以增加计算量为代价的. 从上述的迭代公式可以看出，每迭代一次都要重新计算f x, y 的值，而且迭代又要进行假设干次，计算相当的复杂 . 为此，有没有比拟简便的计算方法呢？下面给出改良的欧拉方法.1.2.3 改良的欧拉方法由前面的讨论可知， 梯形法计算相对复杂， 现对上面的梯形法进行简化， 具体方法是只计算一两次就转入下一步的计算，先用欧拉公式〔 3〕求得一个初步的近似解 y n 1 ，称为预测值，再利用公式〔 6〕把它校正一次，这样建立的预测 - 校正系统通常称为改良的欧拉公式 . 具体公式如下yn 1y nh f x n , y n f x n 1, y n hf x n , y n2〔8〕改良的欧拉法与梯形法一样，是二阶方法.1.2.4 Runge-Kutta方法由前面讨论可知，从〔 4〕式可以看出，假设要使公式阶数提高， 就必须使右端积分的数值求积公式精度提高， 它必然要增加求积积点，为此将〔 4〕式的右端用求积公式表示为xn 1r〔9〕f x, y x dx hc i f x ni h, y x ni h ，x ni1一般来说，点数 r 越多，精度越高，上式右端相当于增量函数 x, y, h ，为得到便于计算的显式方法，将公式〔 9〕表示为：yn 1y n h x n , y n , h ,〔10〕 其中rx n , y n , hc i K ii 1K 1f x n , y n〔11〕i 1K if x n i h, y nh ijKj, i 2, rj 1这里 c i , i , ij 均为常数 . c i 为加权因子， K i 为第 i 段斜率，共有 r 段. 我们把〔 10〕和 〔11〕称为 r 级显式 Runge-Kutta 法，简称为 R-K 方法 . 下面给出其中最经典最常用的一个公式：yynh K 1 2K 2 2K 3 K 4,n 16K 1 f x n , y n ,K 2f x nh, y nhK 122K 3 f x nh, y nhK 2 ,2 2K 4 f x nh, y n hK 3 .〔12〕Runge-Kutta 方法作为一种重要的单步方法，具有很高的实用价值，它关于初值是稳定的，其解连续地依赖于初值，是一类便于应用的单步法，为了计算y n 1 ，只用到前面一步的值 y n 即可，因此每步的步长可以独立取定. 常用的 Runge-Kutta方法精度较高， 为了到达预定的精度， 与欧拉方法与梯形法相比， 步长 h 可取得大些，求解 区间上的总步数可以少 些 . 但 Runge-Kutta 方法也有些 缺点 ，比方 四阶 Runge-Kutta 方法每算一步需要四次计算 f x, y 的值，计算量较大〔对于复杂的f x, y 而言〕．2 数值方法的应用实例 [5-9]y10 y 例 1对于初值问题，分别用欧拉法、 改良的欧拉法， 梯形法求y 1 的y 01近似值 .解：易得该方程的解析解y xe 10x ， y 1，为比拟，将按不同数值计算方法所得结果列表如下：表 1 三种不同方法的数值结果欧拉法 改良的欧拉法 梯形法-1 112.6561 E -005 4.6223 E -005 4.5026 E -0054.3717 E -005 4.5408 E -005 4.5396 E -0054.5173 E -005 4.5400 E -005 4.5400 E -005图 1三种不同方法数值解与精确解的误差曲线从表 1 中可以看出：当 h0.2 时，三种方法均不稳定，计算结果严重偏离精确值；h时，改良后的欧拉和梯形法均稳定，但欧拉法效果很差；当时，三种方法均稳定，但精确度有区别．可以看出， h 越小，计算结果越好，要想计算结果充分接近于解析解还须取较小的 h 值．图 1 反映的步长 h 0.01 时，三种数值方法的所得数值解与解析解在[0,1] 区间的误差曲线，由图可知，在步长相同的情况下，梯形法的精确度略高于改良的欧拉法；改良的欧拉法和梯形法精确度都明显高于欧拉法．例 2 用欧拉法、改良的欧拉法和 Runge-Kutta 法求解初值问题并比拟三种方法的结果．解：方程为 n 1 的伯努利方程，可求得解析解为现用 MATLAB软件编程，用题目要求的方法求解，可得如下列图示结果：图 2 〔a〕步长为 0.2 时 R-K 法和解析解比拟图 2 〔b〕步长为 0.2 时改良的 Euler 法和解析解比拟图 2 〔c〕步长为 0.2 时欧拉法和解析解比拟上图 2(a) ，(b) ，(c) 描述的是步长为0.2 时，用欧拉法、改良的欧拉法， Runge-Kutta 法求解方程所得的数值解与解析解之间的比照图．由图可知，Runge-Kutta 法所得数值解曲线和解析解曲线吻合的很好，改良的欧拉法和欧拉法随着计算的进行，数值解和解析解之间误差逐步增大，但改良的欧拉法效果要好于欧拉法．图 3 (a)步长为时Euler法和解析解比较图 3 (b)步长为时改良的Euler法和解析解比拟图 3 (c) 步长为 0.1 时 Runge-Kutta 法和解析解比拟上图 3 (a) ，(b) ，(c) 描述的是步长为0.1 时，用欧拉法、改良的欧拉法， Runge-Kutta 法求解方程所得的数值解与解析解之间的比照图．由图可知，改良的欧拉法和Runge-Kutta 法所得数值解曲线和解析解曲线吻合的很好，而欧拉法随着计算的进行数值解和解析解之间误差逐步增大．相应的程序如下：主程序x=0:0.2:1;jxj=exp(2*x).*(1./exp(4*x) + (2*x)./exp(4*x)).^(1/2);y=Euler(@ff,0,1,0.2,1);gy=geuler(@ff,0,1,0.2,1);Ry=RK(@ff,0,1,0.2,1);figure(1);plot(x,jxj,x,Ry,'*');figure(2);plot(x,jxj,x,gy,'*');figure(3);plot(x,jxj,x,y,'*' )欧拉法程序function y=Euler(f,a,b,h,y0)n=(b-a)/h;x=a:h:b;y=zeros(n+1,1);y(1)=y0;for i=1:ny(i+1)=y(i)+h*feval(f,x(i),y(i));end改良的欧拉法程序function gy=geuler(f,a,b,h,y0)n=(b-a)/h;x=a:h:b;y=zeros(n+1,1);y(1)=y0;for i=1:nyp=y(i)+h*feval(f,x(i),y(i));yc=y(i)+h*feval(f,x(i+1),yp);y(i+1)=(yp+yc)/2;endgy=y;Runge-Kutta 法程序function Ry=RK(f,a,b,h,y0)n=(b-a)/h;x=a:h:b;y=zeros(n+1,1);y(1)=y0;for i=1:nk1=feval(f,x(i),y(i));k2=feval(f,x(i)+h/2,y(i)+h*k1/2);k3=feval(f,x(i)+h/2,y(i)+h*k2/2);k4=feval(f,x(i+1),y(i)+h*k3);y(i+1)=y(i)+h*(k1+2*k2+2*k3+k4)/6;endRy=y;3微分方程数值解法在实际生活中的应用3.1 应用实例 : 耐用消费新产品的销售规律模型一种新产品进入市场以后，常常会经历销售量首先慢慢增加然后逐渐慢慢下降的一个过程，由此给出的时间—销售坐标系下的曲线称为产品的生命曲线 , 它的形状是钟形的 . 不过对于较耐用的消费品，情况会有所不一样，它的生命曲线会在开始有个小小的顶峰，之后是段比拟平坦的曲线，先下降，之后再上升，然后到达顶峰，因此它是双峰型的曲线 .如何理解这种和传统产品的生命曲线理论相冲突的现象？澳大利亚学者斯蒂芬斯与莫赛观察到的购置耐用消费产品的人大概可分为两类：其一是非常善于接受新的事物，称其为“创新型〞消费者，他们会经常从产品广告，制造商给出产品的说明书与商店样品中了解了产品功能与性能之后，再决定其否购置；其二是消费者相对保守，称其为“模仿型〞消费者，他们往往会根据大局部已经购置产品的消费者实际使用的经验而提供的信息来决定是否购置产品，下面的例子经过分析，建立相应的数学模型，对这种现象给出了科学解释.模型的建立将顾客获得信息大致可分成两类，其一称之为“搜集型〞，源自于产告说明、广告，“创新型〞顾客获得这类消息后就可作出其是否购置；第二类称之为“体验型〞，即消费者使用之后会获得真实体验，常常以口头的形式散播，“模仿型〞类顾客会在获得这种信息之后才能决定是否购置 .设 K 是潜在用户的总数，K1与 K 2分别为“创新型〞与“模仿型〞的人数. 再设 N t 是时刻t已经购置的商品顾客的人数，而N 1 t 与 N 2 t 分别为“创新型〞与“模仿型〞的顾客的人数, 再设 A1 t 是时刻t中已获得“搜集型〞信息的人数，由于该局部的信息可直接由外部获得，或者可从已获得该类信息的人群中获得，因此类似巴斯模型，从而建立如下方程：dA1t12 A1t , A 0 0,1,20，dt K1 A1 t获得“搜集型〞信息的“创新型〞消费者决定其是否购置的行为，满足如下方程：对于“模仿型〞的顾客，可从已经购置该产品“创新型〞或者“模仿型〞的顾客中获取信息，于是有在这里，忽略消费者购置产品后需一段短暂的使用后才会散播体验信息的滞后作用．综上，斯蒂芬斯—莫赛模型是一常微分方程组的初值问题：而 N t N1 t N2 t 为时刻t购置该商品的总人数 .模型的求解对于斯蒂芬斯—莫赛模型中N2 t的解析解那么不能求出，于是可以用四阶Runge Kutta 公式求得，且在它的精度要求到达很高情形下求出N2 t .用 MATLAB软件求解，相应的程序及结果如下function RK=RKFC(fc,a,b,h,y0)n=(b-a)/h;x=a:h:b;m=length(y0);y=zeros(n+1,m);y(1,:)=y0;for i=1:nk1=feval(fc,x(i),y(i,:));k2=feval(fc,x(i)+h/2,y(i,:)+h*k1/2);k3=feval(fc,x(i)+h/2,y(i,:)+h*k2/2);k4=feval(fc,x(i+1),y(i,:)+h*k3);y(i+1,:)=y(i,:)+h*(k1+2*k2+2*k3+k4)/6;endRK=y;function f=FC(x,y)k1=50; k2=70;al=0.0013;be=0.0013;ga=0.0015;f=[(k1-y(1))*(al+be*y(1)),ga*(k2-y(2))*(y(1)+y(2))];x=0:0.3:24;微分方程数值解RK=RKFC(@FC,0,24,0.3,[0,0])figure(1) ；plot(x,RK(:,1),'+' ,x,RK(:,2), '*' ) ；legend( 'N1(t)' , 'N2(t)',2) figure(2) ；plot(x,RK(:,1)+RK(:,2),'+' ) ；legend( 'N(t)',2)图 4N1 t , N2 t 与时间关系图图 5N t 与[0,25]时间段关系图由此例可以看出，微分方程数值解在实际生活有着广泛的应用，而数值解法中的Runge-Kutta 方法是一种很重要且应用很广泛的算法 .结语微分方程数值解是求解微分方程的一种很重要且应用范围很广的方法，而常用的数值解法如欧拉法、改良的欧拉法、梯形法和Runge-Kutta 方法在一定程度上都有自己的优缺点，理解和掌握各种方法的应用范围，用 MATLAB编制各种方法求解实际问题的通用程序，对用微分方程数值解理论解决现实生活中的实际问题有很重要的意义．参考文献[1]李庆扬 , 王能超 , 易大义 . 数值分析 [M]. 北京 : 清华大学出版社 , 清华大学出版社 ,2001.[2]冯康 . 数值计算方法 [M]. 杭州 : 浙江大学出版社 ,2003.[3]封建湖 . 计算方法典型题分析解集 [M]. 西安 : 西北工业大学出版社 ,2000.[4]胡建伟 , 汤怀民 . 微分方程数值解法 ( 第二版 )[M]. 北京 : 科学出版社 ,2007.[5]王能超 . 计算方法简明教程 [M]. 北京 : 高等教育出版社 ,2004.[6]Nash S G.A history of scientific computing.[M] New York:ACM Press,1990.[7]于丽妮 . ODE 问题解析解及数值解的 matlab 实现 [J]. 电脑知识与技术 .2021,8(14):3303-3305.[8]霍晓成 . 常微分方程数值解法的研究 [J]. 临沂师范学院学报 . 2021,(6):19-23.[6]王国立 ,陈瑛.非线性微分方程迭代算法及其在物理学中的应用[J].长春师范学院学报(自然科学版).2006, 25(2):10-12.。

微分方程数值解法实验报告姓名： 班级： 学号：一：问题描述求解边值问题：()2(sin cos cos sin (0,1)(0,1)0,(,)x y u e x y x y G u x y G ππππππ+⎧⎫∆=+⎪⎪∈=⨯⎨⎬⎪⎪=∈∂⎩⎭（x,y) 其精确解为)sin()sin(),()(y x e y x u y x πππ+=问题一：取步长h=k=1/64,1/128,作五点差分格式，用Jacobi 迭代法，Gauss_Seidel 迭代法，SOR 迭代法（w=1.45）。

求解差分方程，以前后两次重合到小数点后四位的迭代值作为解的近似值，比较三种解法的迭代次数以及差分解)128/1,64/1)(,(=h y x u h 与精确解的精度。

问题二：取步长h=k=1/64,1/128,作五点差分格式，用单参数和双参数PR 法解差分方程，近似到小数点后四位。

与SOR 法比较精度和迭代步数。

问题三：取步长h=k=1/64,1/128,作五点差分格式，用共轭梯度法和预处理共轭梯度法解差分方程，近似到小数点后四位。

与SOR法与PR 法比较精度和迭代步数。

二．实验目的：分别使用五点差分法（Jacobi 迭代，Gauss_Seidel 迭代，SOR迭代），PR 交替隐式差分法（单参数，双参数），共轭梯度法，预共轭梯度法分别求椭圆方程的数值解。

三．实验原理：(1) Jacobi 迭代法设线性方程组(1)的系数矩阵A 可逆且主对角元素均不为零,令 并将A 分解成(2) 从而(1)可写成令其中. (3) 以为迭代矩阵的迭代法(公式)(4)称为雅可比(Jacobi)迭代法(公式),用向量的分量来表示,(4)为(5) 其中为初始向量. (2) Guass-Seidel 迭代法由雅可比迭代公式可知,在迭代的每一步计算过程中是用的全部分量来计算的所有分量,显然在计算第i 个分量时,已经b Ax =nn a ,...,a ,a 2211()nn a ,...,a ,a diag D 2211=()D D A A +-=()b x A D Dx +-=11f x B x +=b D f ,A D I B 1111--=-=1B ()()111f x B x k k +=+⎩⎨⎧[],...,,k ,n ,...,i x a b a x n ij j )k (j j i i ii )k (i 21021111==∑-=≠=+()()()()()Tn x ,...x ,x x 002010=()k x ()1+k x ()1+k i x计算出的最新分量没有被利用，从直观上看,最新计算出的分量可能比旧的分量要好些.因此，对这些最新计算出来的第次近似的分量加以利用,就得到所谓解方程组的高斯—塞德（Gauss-Seidel ）迭代法.把矩阵A 分解成(6)其中,分别为的主对角元除外的下三角和上三角部分,于是,方程组(1)便可以写成 即其中(7)以为迭代矩阵构成的迭代法(公式)(8)称为高斯—塞德尔迭代法(公式),用 量表示的形式为(3) SOR 迭代(4) 交替方向迭代法（PR 法）迭代格式为：()()1111+-+k i k x ,...,x 1+k ()1+k x ()1+k j x U L D A --=()nn a ,...,a ,a diag D 2211=U ,L --A ()b Ux x L D +=-22f x B x +=()()b L D f ,U L D B 1212---=-=2B ()()221f x B x k k +=+⎩⎨⎧[],...,,k ,n ,,i x a x a b a x i j n i j )k (j ij )k (j ij i ii )k (i 21021111111==∑∑--=-=+=++Λ))1(()(1D R L D T ωωω-+-=-b )(1--=L D d ωωhu πμωcos )11/(22opt =-+=2121,,1,1,1,,122L L L L u u u L u u u j i j i j i j i j i j i +==+-=+-+-+-对于单参数PR 法，对于多参数，(5) 共轭梯度法 算法步骤如下： ［预置步］任意，计算，并令取：指定算法终止常数，置，进入主步；［主步］ （１）如果，终止算法，输出；否则下行；（２）计算：（３）计算：（４）置，转入（１）．(6) 预共轭梯度法b uL I uL I b u L I uL I k k k k k k k k k k ττττττ+-=++-=++++211122211)()()()(hh optπτsin 22=2sin a ....2,1)11(421k 221h k a h k πρρτ==+-=--其中［预置步］任意，计算，并令取：指定算法终止常数，置，进入主步；［主步］（１）计算：，（２）如果，转入（3）．否则，终止算法，输出计算结果（３）计算：（４）置，转入（1）注：在算法［主步］中，引入变量，及，可以简化计算。

微分方程的常用数值解法摘要：微分方程是数学中的一种重要的方程类型，它能描述自然现象和工程问题中的许多变化规律。

但是大多数微分方程解法是无法用解析的方式求解的，因此需要借助数值解法来近似求解。

本文将介绍微分方程的常用数值解法。

关键词：欧拉方法；龙格-库塔方法；微分方程；常用数值解法一、微分方程数值解方法微分方程数值解法是数学中的重要部分。

欧拉方法、龙格-库塔方法和二阶龙格-库塔方法是常用的微分方程数值解法，下面就分别介绍这三种方法。

（一）欧拉方法欧拉方法是解初值问题的一种简单方法，它是欧拉用的第一种数值方法，也叫向前欧拉法。

欧拉方法是利用微分方程的定义式y’=f(x, y)，将它带入微分方程初值问题y（x_0）=y_0中，以y_0为初始解，在每一步上通过沿着切线的方法进行估计并推进新的解y_{i+1}：y_i+1=y_i+hf(x_i,y_i)其中，x_i和y_i是我们知道的初始条件，h是求解过程中的步长，f是微分方程右端项。

它是一种时间迭代的算法，易于实现，但存在着精度不高的缺点。

（二）龙格-库塔方法龙格-库塔方法是一种经典迭代方法，也是近代微分方程数值解法发展的里程碑之一。

龙格-库塔方法的主要思想是利用规定的阶码及阶向量，通过递推求解微分方程数值解的近似值。

龙格-库塔方法的方式不同，其步骤如下：第一步：根据微分方程，计算出在x_i和y_i的值。

第二步：在x_i处对斜率进行估计，并利用这个斜率来求解下一步所需的y_i+1值。

第三步：使用x_i和y_i+1的值来重新估计斜率。

第四步：使用这个新的斜率来更新y_i+1的值。

（三）二阶龙格-库塔方法二阶龙格-库塔方法是龙格-库塔方法的一种变体，它根据龙格-库塔方法的思想，使用更好的步长来提高数值解的精度。

二阶龙格-库塔方法的基本思路是，在第一次迭代时使用一个阶段小一半的y_i+1，然后使用这个估算值来计算接下来的斜率。

通过这种方法，可以提高解的精度。

二阶龙格-库塔方法的步骤如下：第一步：计算出初始阶段的y_i+1值。

随机微分方程(stochastic differential equation,sde) 1. 引言1.1 概述随机微分方程（Stochastic Differential Equation，SDE）是一类描述随机现象的微分方程。

相比于传统的确定性微分方程，SDE中包含了一个或多个随机项，能够更准确地描述现实世界中的不确定性和变动性。

SDE在各个领域中广泛应用，特别是金融学、物理学和生物学等领域。

1.2 文章结构本文将从以下几个方面介绍随机微分方程及其应用：定义与基本概念、解随机微分方程的方法与技巧，以及在实际问题中的应用。

具体可以分为三个主要部分：引言、主体内容和结论展望。

1.3 目的本文旨在介绍随机微分方程的基本概念、解法和应用，并探讨其在金融学、物理学和生物学等领域中的实际应用。

通过对随机微分方程的深入了解，读者可以更好地理解和利用该方法来解决实际问题，并对未来研究提出展望。

以上为“1. 引言”部分的内容。

2. 随机微分方程的定义与基本概念2.1 随机过程简介随机过程是一类描述随着时间推移而随机变化的数学模型。

它可以看作是时间参数上的一族随机变量的集合。

随机过程常用于描述具有随机性质的现象，如金融市场中的股票价格、天气预报中的温度变化等。

2.2 随机微分方程的定义随机微分方程是一类描述含有随机项（通常为噪声）的微分方程。

它通常采用以下形式表示：dX(t) = a(X(t), t)dt + b(X(t), t)dW(t)其中，X(t)是未知函数，a(X(t), t)和b(X(t), t)是已知函数，dW(t)表示Wiener 过程（也称为布朗运动或白噪声）。

这个方程表示了X在无穷小时间段dt内发生微小变化dX(t)，其中包含一个确定性项a(X(t), t)dt和一个随机项b(X(t), t)dW(t)。

2.3 常见的随机微分方程模型在实际应用中，有许多不同类型的随机微分方程模型被广泛使用。

- Ornstein-Uhlenbeck 过程：该模型描述了维持平衡状态的粒子在受到随机扰动时的演化过程。