随机微分方程数值解法.
微分数值解法
微分数值解法分别用欧拉方法、改进的欧拉方法、二阶龙格-库塔方法、三阶龙格-库塔方法、四阶龙格-库塔方法,计算下列初值问题,并与精确解对比。
微分方程的初值问题:⎪⎩⎪⎨⎧==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)四阶龙格-库塔方法:。
微分方程数值解第一章答案
u(t) hf (t, u(t))
23
数值方法的基本问题
• 截断误差(局部、整体) • 相容性 • 收敛性 • 稳定性
24
局部截断误差
设u(t)是初值问题(1)的解, 在[t,t+h]上定义算子
R(t, u; h) u(t h) u(t) hf (t, u(t))
通解 — 解中所含独立的任意常数的个数与方程 的阶数相同. (微分方程的绝大部分解)
特解 — 不含任意常数的解.
例y Cex是方程y y的通解.
例y 2ex是方程y y的特解.
11
➢ 定解条件: 初值问题和边值问题
定解条件 — 确定通解中任意常数的条件.
1) n 阶方程的初始条件(或初值条件):
那么, n称为整体截断误差
与局部截断误差不同, 此时
u(tk ) uk,k 0,1,L n 未必成立, 且一般 u(tk ) uk , k 0,1,L n
26
截断误差
• 局部截断误差Rn:假设第n步精确计算的前 提下,计算解un+1和精确解u(tn+1)的误差
• 整体截断误差n:在考虑误差累积的效应
北京·中国地质大学
China University of Geosciences,Beijing
《微分方程数值解法》
教材: 微分方程数值方法
(第二版), 胡健伟,汤怀民著, 科学出版社, 2007,2
参考书: 微分方程数值解法 李荣华等编, 高教出版社 1
参考书:
微分方程数值解法 李荣华等编, 高教出版社
u' f (t, u),
微分方程数值解法
微分方程数值解法微分方程数值解法微分方程数值解法【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方法精度高。
微分方程和偏微分方程的数值解法
描述金融衍生品的定价过程,如布莱克-舒尔斯模型就是一个偏微分方程。通过求解该 方程,可以得到期权的理论价格以及相应的风险参数。
投资组合优化
在投资组合理论中,常使用微分方程来描述资产价格的动态变化和投资者的风险偏好。 通过求解这些方程,可以得到最优的投资组合配置策略以实现风险与收益的平衡。
数值解法需要保证稳定性和收敛 性,即当离散间隔趋近于零时, 数值解应趋近于真实解。
02
常微分方程的数值解法
欧拉方法
基本思想
通过逐步逼近的方式,利用已知点的信 息来推测下一个点的信息。
公式推导
基于泰勒级数展开,忽略高阶项得到近 似公式。
优缺点
简单易懂,但精度较低,仅适用于简单 问题。
改进方法
采用改进的欧拉方法或预估-校正法提 高精度。
物理问题中的微分方程和偏微分方程
牛顿第二定律
描述物体运动的基本定律,可以 表示为二阶常微分方程。通过求 解该方程,可以得到物体的位移 、速度和加速度等运动学量。
热传导方程
描述热量在物体内部传递的过程 ,是一个偏微分方程。通过求解 该方程,可以得到物体内部的温 度分布以及热量的传递速率。
波动方程
描述波动现象(如声波、光波等 )的传播过程,是一个二阶偏微 分方程。通过求解该方程,可以 得到波的传播速度、振幅、频率 等波动特性。
工程问题中的微分方程和偏微分方程
结构力学中的弹性力学方程
描述结构在受力作用下的变形和应力分布,是一个偏微分方程。通过求解该方程,可以得到结构的位移、应 力和应变等力学量,为工程设计提供重要依据。
流体力学中的纳维-斯托克斯方程
描述流体运动的基本方程,是一个偏微分方程。通过求解该方程,可以得到流体的速度、压力和温度等流场 特性,为流体机械设计和优化提供指导。
微分方程数值解法答案
微分⽅程数值解法答案包括基本概念,差分格式的构造、截断误差和稳定性,这些内容是贯穿整个教材的主线。
解答问题关键在过程,能够显⽰出你已经掌握了书上的内容,知道了解题⽅法。
这次考试题⽬的类型: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.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) 共轭梯度法 算法步骤如下: [预置步]任意,计算,并令取:指定算法终止常数,置,进入主步;[主步] (1)如果,终止算法,输出;否则下行;(2)计算:(3)计算:(4)置,转入(1).(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 πρρτ==+-=--其中[预置步]任意,计算,并令取:指定算法终止常数,置,进入主步;[主步](1)计算:,(2)如果,转入(3).否则,终止算法,输出计算结果(3)计算:(4)置,转入(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)
随机微分方程(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 过程:该模型描述了维持平衡状态的粒子在受到随机扰动时的演化过程。
随机微分方程数值解法
方程(9)即为Stratonovich型随机微分方程。 注:1)Itó型随机微分方程(6)和Stratonovich型随机微分方程(9) 是可以相互转换的。在标量情形下,对方程(6)令
plot([0:dt:T],[0,W],’r-’) %绘图 xlabel(’t’,’FontSize’,16) ylabel(’W(t)’,’FontSize’,16,’Rotation’,0)
1.2 随机积分
随机积分分为Itó型随机积分和Stratonovich型随机积分。以
下假设Wiener过程 W(t),t0定义在概率空间 (,F,P)上,
0 t 0 t 1 t 2 t n t ,
令 t k t k t k 1 ( 1 k n ) , m 1 k a x n t k ,
若随机变量序列
n
X (tk 1 )(W (tk) W (tk 1 )),n 1 ,2 ,3
(4)
k 1
均方收敛于唯一极限,则称
f ( t , x ) f ( t , y ) g ( t , x ) g ( t , y ) L 2 x y , x R , 且有E y0 2 , 则方程 (6)存在唯一解且E y(t)2 。
定义 2.1 (强收敛性) 若存在常数 C 0(与 h 独立), 0 ,使得
E (y ( tn ) y n ) C h p ,h ( 0 ,) ,
设 是正整数,
利用随机
Taylor展开式和Itó公式,可以得到:
y ( t n 1 ) y ( t n ) I 0 f ( y ( t n ) ) I 1 g ( y ( t n ) ) I 1 1 L 1 g ( y ( t n ) ) I 0 0 L 0 f ( y ( t n ) ) R ,( 1 1 ) 其中R 是余项,算子 L 0 和 L 1 分别为
微分方程数值解第一章答案
微分方程数值解法 李荣华等编, 高教出版社
• 课堂授课+计算实验 • 考核方式: 平时作业+课堂+期末考试 • 任课教师 •
1
教学内容
• 第一章、常微分方程的数值解法 • 第二章、椭圆型方程的差分方法 • 第七章、椭圆型方程的有限元方法 • 第四章、抛物型方程的差分方法 • 第五章、双曲型方程的差分格式
x (a bx)t x ' ax bx2 x
Logistic方程
14
常微分方程举例3
问题1.3 并不是所有的方程可以用初等积分法求出其 解, 例如形式上很简单的里卡蒂(Riccati)方程
x' t2 x2
不能用初等函数表示通解. 寻求方程非解析函数的其它形式解, 显得非常必要。 而数值求解就是其重要的一个方法
问题得真解;即收敛性问题 ② 误差估计 ③ 产生得舍入误差,在以后得各步计算中,是否会
无限制扩大;稳定性问题
32
数值求解微分方程过程示意
区域剖分
微分方程离散
微
初始和边界条件处理
分
解的存在性、唯一性
方
离散系统的 性态研究
解的收敛性和收敛速度
程
解的稳定性
递推计算或解线 性代数方程组
得到数值解
33
作业
18
举例2
• P55 习题1 利用Euler方法求数值解 初值问题u' 1 u, u(0) 1 2 步长h=0.1, 解区间[0,1]
• 绘制折线,与真解比较
19
Matlab实现 u=null(1);h=0.1;u0=1; u(1)=u0+h*0.5*u0; for n=1:9
u(n+1)=u(n)+h*0.5*u(n); end t=0:0.1:1;un=[u0,u]; plot(t,un,'ro','Linewidth',2) ut=exp(0.5*t); hold on plot(t,ut,'Linewidth',2)
微分方程数值解实验报告
微分方程数值解法课程设计报告班级:姓名:学号:成绩:2017年 6月 21 日摘要自然界与工程技术中的很多现象,可以归结为微分方程定解问题。
其中,常微分方程求解是微分方程的重要基础内容。
但是,对于许多的微分方程,往往很难得到甚至不存在精确的解析表达式,这时候,数值解提供了一个很好的解决思路。
,针对于此,本文对常微分方程数值解法进行了简单研究,主要讨论了一些常用的数值解法,如欧拉法、改进的欧拉法、Runge—Kutta方法、Adams法以及椭圆型方程、抛物型方程的有限差分方法等,通过具体的算例,结合MATLAB求解画图,初步给出了一般常微分方程数值解法的求解过程。
同时,通过对各种方法的误差分析,让大家对各种方法的特点和适用范围有一个直观的感受。
关键词:微分方程数值解、MATLAB目录摘要 (2)目录 (3)第一章常微分方程数值解法的基本思想与原理 (4)1.1常微分方程数值解法的基本思路 (4)1.2用matlab编写源程序 (4)1.3常微分方程数值解法应用举例及结果 (5)第二章常系数扩散方程的经典差分格式的基本思想与原理 (6)2.1常系数扩散方程的经典差分格式的基本思路 (6)2.2 用matlab编写源程序 (7)2.3常系数扩散方程的经典差分格式的应用举例及结果 (8)第三章椭圆型方程的五点差分格式的基本思想与原理 (10)3.1椭圆型方程的五点差分格式的基本思路 (10)3.2 用matlab编写源程序 (10)3.3椭圆型方程的五点差分格式的应用举例及结果 (12)第四章总结 (12)参考文献 (12)第一章常微分方程数值解法的基本思想与原理1.1常微分方程数值解法的基本思路常微分方程数值解法(numerical methods forordinary differential equations)计算数学的一个分支.是解常微分方程各类定解问题的数值方法.现有的解析方法只能用于求解一些特殊类型的定解问题,实用上许多很有价值的常微分方程的解不能用初等函数来表示,常常需要求其数值解.所谓数值解,是指在求解区间内一系列离散点处给出真解的近似值.这就促成了数值方法的产生与发展.1.2用matlab编写源程序龙格库塔法:M文件:function dx=Lorenz(t,x)%r=28,sigma=10,b=8/3dx=[-10*(x(1)-x(2));-x(1)*x(3)+28*x(1)-x(2);x(1)*x(2)-8*x(3)/3];运行程序:x0=[1,1,1];[t,y]=ode45('Lorenz',[0,100],x0);subplot(2,1,1) %两行一列的图第一个plot(t,y(:,3))xlabel('time');ylabel('z');%画z-t图像subplot(2,2,3) %两行两列的图第三个plot(y(:,1),y(:,2))xlabel('x');ylabel('y'); %画x-y图像subplot(2,2,4)plot3(y(:,1),y(:,2),y(:,3))xlabel('x');ylabel('y');zlabel('z');%画xyz图像欧拉法:h=0.010;a=16;b=4;c=49.52;x=5;y=10;z=10;Y=[];for i=1:800x1=x+h*a*(y-x);y1=y+h*(c*x-x*z-y);z1=z+h*(x*y-b*z);x=x1;y=y1;z=z1;Y(i,:)=[x y z];endplot3(Y(:,1),Y(:,2),Y(:,3));1.3常微分方程数值解法的应用举例及结果应用举例:⎪⎪⎪⎩⎪⎪⎪⎨⎧-=--=--=)()()()()()()()()())()(()(t bz t y t x dt t dz t z t x t y t rx dt t dy t y t x a dt t dx a=10,b=8/3,0<r<+∞,当1<r<24.74时,Lorenz 方程有两个稳定的不动点c()1(-r b ,)1(-r b ,r-1)和c '(-)1(-r b ,-)1(-r b ,r-1),一个稳定的不动点0=(0,0,0),当r>24.74时,c 和c '都变成不稳定的,此时存在混沌和奇怪吸引子。
微分方程的基本概念和解法技巧
微分方程的基本概念和解法技巧微分方程是数学中重要的一种方程,它涉及到函数与它的导数之间的关系。
在物理学、工程学、经济学等领域中,微分方程广泛应用于描述各种变化和运动的规律。
了解微分方程的基本概念和解法技巧,对于理解和解决实际问题具有重要意义。
本文将介绍微分方程的基本概念以及一些常见的解法技巧。
一、微分方程的基本概念1. 定义:微分方程是含有未知函数及其导数的方程。
一般形式可以表示为 F(x, y, y', y'', ...) = 0,其中 y 是未知函数。
2. 阶数:微分方程的阶数是指该方程中导数的最高阶数。
常见的阶数有一阶、二阶和高阶微分方程。
3. 解:微分方程的解是满足方程的函数。
一般来说,一个微分方程可以有无穷多个解。
4. 初值问题:初值问题是求解微分方程时给定一个或多个初始条件,根据这些条件确定方程的解。
初值问题通常涉及到一个点上的初始状态。
5. 常微分方程和偏微分方程:常微分方程只涉及到一个自变量,而偏微分方程则涉及到多个自变量。
常微分方程的解是一类函数,而偏微分方程的解是一个函数族。
二、微分方程的解法技巧1. 变量可分离法:适用于可以将微分方程的变量分离开的情况。
通过将方程两边同时乘以不同变量的函数,使得方程可以变为两个积分的形式,从而得到解。
2. 齐次方程法:适用于可以通过变量代换将微分方程化为齐次方程的情况。
齐次方程中的未知函数可以表示为一个比值函数,通过变量代换后,方程可以化为一个仅依赖于一个变量的方程,从而得到解。
3. 一阶线性常微分方程:适用于形如 y' + p(x)y = q(x) 的一阶线性常微分方程。
通过乘以一个适当的积分因子将方程化为可积形式,然后求解积分得到方程的解。
4. 常系数线性微分方程:适用于形如 y⁽ⁿ⁾ + aₙy⁽ⁿ⁻¹⁾ + ... + a₁y' + a₀y =g(x) 的常系数线性微分方程。
通过猜测形式,得到特解和齐次方程的通解,从而得到方程的通解。
微分方程的数值解法
4 微分方程的数值解法 4.2 常微分方程边值问题的数值解
方程离散:
Gx 和 Rx 在各节点的值表示为 将各函数 Px 、 Pi Pa ih 、 Gi Ga ih 和 Ri Ra ih 。原方程中各导 数用差商代替有
dy பைடு நூலகம் x h y x h dx 2h
求解思路:区间离散、方程离散 区间离散: 将a到b整个范围内分成n个等距区间,令 h b a / n , 则第i个区间的终点为 xi a ihi 0,1,2,, n ,在该点的y 值可表示为 yi ya ihi 0,1,2,, n 。 求解代数方程组
式中yj,i为第j个变量yj(x)在节点xi处的近似解; n为因变量和方程的个数。
4 微分方程的数值解法 4.2 常微分方程边值问题的数值解 ——二阶常微分方程的有限差分法
二阶常微分方程的边值问题: d 2 y dy 2 P x G x y R x dx dx ya A, yb B
4 微分方程的数值解
4.1.4 常微分方程组初值问题的数值解 一个自变量,m个因变量组成的一个常微分方程组
f1 x, y1 , y2 , , ym y1 y f x, y , y , , y 2 2 1 2 m f m x, y1 , y2 , , ym ym y x y 2 0 20 y m x0 y m 0
4 微分方程的数值解法 基本概念
微分方程的初值问题: 求解微分方程时,必须有一些已知条件。若所给 的已知条件为某特定点上各阶因变量的值,此类问题 为初值问题。
dCA kCA dt t 0, C A C A0
微分方程数值解法(戴嘉尊)习题解答
n+1
n
n
n
n
n
n
n
(n=0,1,2,…………9)
y y y x y x = + h [− + +1 + (− + +1)]
n+1
n2
n
n
n+1
n+1
y x y x h
=(1- )
2
+ h ( + 1) + h (-
n2 n
2
+
n+1
n+1 + 1 )
y x y x =0.95 + 0.05( +1) + 0.05(− + +1)
h 2
f
(xn , y(xn )) −
yn
−
h 2
f (xn , yn )) |
≤
Lh[|
y(xn ) −
yn
|+
h 2
|
f (xn , y(xn )) −
h 2
f (xn , yn )) |]
≤
Lh(|
εn
|
+
h 2
L
|
y ( xn
)
−
yn
|]
≤
Lh(1 +
h 2
L)
|
εn
|
综上有:
|
ε n+1
|≤|
(xn ,
yn ))]dx
∫=
xn+1 [y′(x) −
xn
y′(x +
h 2
)]
+
数值计算方法实验报告
数值计算方法实验报告一、实验目的本实验旨在通过数值计算方法的实验操作,深入理解数值计算方法的原理与应用,掌握数值计算方法的相关技能,提高数值计算方法的实际应用能力。
二、实验内容1.数值微积分2.数值代数3.数值微分方程4.数值线性代数5.数值优化6.数值统计分析7.数值随机模拟8.数值傅立叶分析9.数值偏微分方程三、实验步骤1.数值微积分:通过不同的数值积分方法,计算给定函数的定积分值,并对不同数值积分方法的误差进行分析。
2.数值代数:通过使用线性代数方法,求解给定的线性方程组,并分析不同线性方程组求解方法的优劣。
3.数值微分方程:通过使用常微分方程数值解法,求解给定的微分方程,并比较不同求解方法的精度和稳定性。
4.数值线性代数:通过使用特征值分解方法,对给定的矩阵进行特征值分解,并分析不同特征值分解方法的优缺点。
5.数值优化:通过使用不同的优化方法,求解给定的优化问题,并比较不同的优化方法的效率和精度。
6.数值统计分析:通过使用不同的统计分析方法,对给定的数据进行统计分析,并分析不同的统计方法的优缺点。
7.数值随机模拟:通过使用随机模拟方法,模拟给定的概率分布,并分析不同随机模拟方法的效率和精度。
8.数值傅立叶分析:通过使用傅立叶分析方法,对给定的信号进行频谱分析,并分析不同的傅立叶分析方法的优缺点。
9.数值偏微分方程:通过使用偏微分方程数值解法,求解给定的偏微分方程,并比较不同求解方法的精度和稳定性。
四、实验结果与分析本实验中,通过对不同的数值计算方法的实验操作,我们可以更深入地理解数值计算方法的原理与应用,并掌握数值计算方法的相关技能,提高数值计算方法的实际应用能力。
同时,通过实验结果的分析,我们可以更好地比较不同数值计算方法的优缺点,为实际应用提供参考依据。
五、实验总结本实验旨在通过数值计算方法的实验操作,深入理解数值计算方法的原理与应用,掌握数值计算方法的相关技能,提高数值计算方法的实际应用能力。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
的数值算法。
N
h (T t0 ) / N ,tn t0 nh,
设 是正整数,
利用随机
Taylor展开式和Itó公式,可以得到:
y(tn1) y(tn) I0 f (y(tn)) I1g(y(tn)) I11L1g(y(tn)) I00L0 f (y(tn)) R, (11) 其中R 是余项,算子 L0和 L1 分别为
随机微分方程数值解法
2013年11月18日
随机微分方程数值解法
1.随机微分方程概述
1.1 布朗运动介绍 1.2 随机积分 1.3 两种形式的随机微分方程
2.随机微分方程数值方法介绍
2.1 随机Taylor展开 2.1 Euler方法 2.2 Milstein方法
3. 数值试验
3.1 精度数值试验 3.2 稳定性数值试验
是可以相互转换的。在标量情形下,对方程(6)令
f (t, y(t )) f (t, y(t )) 1 g (t, y(t ))g(t, y(t )), 2 y
在矢量情形下,令
f i (t, y(t ))
fi
(t,
y(t ))
1 2
m j 1
d k 1
gik y j
(t,
y(t ))g jk (t ,
2) 有些随机微分方程的解析解虽然可以求到,但是形式很复 杂,处理起来很不方便。
3) 在实际应用中,实用的方法是在计算机上进行数值求解,
即不直接求出 y(t ) 的解析解,而是在解所存在的区间上,求得一 系列点 xn(n 1, 2, ) 上的近似值。
2.随机微分方程数值方法介绍
目前随机微分方程的数值求解方法有Euler方法、Milstein方 法 、Runge-Kutta方法等。Runge-Kutta方法的复杂程度比Euler 方法和Milstein方法的程度要高。在实际应用中,一般情况下用 Euler方法和Milstein方法来对模型进行数值模拟。由于Itó型随机 微分方程与Stratonovich型随机微分方程是可以相互相互转化的, 以下介绍求解Itó型随机微分方程(6)的Euler方法和Milstein方法。
y(t )),
其中 i 1, , m. 则方程(6)可以转化为Stratonovich性随机微
分
方程如d下y:(t ) f (t, y(t ))dt g(t, y(t )) dW (t ).
注:1) 大部分随机微分方程的解析解是无法获得的,可以求得解 析解的随机微分方程多为线性随机微分方程。
%布朗运动的模拟
randn('state',100)
% 设置随机数发生器的状态
T=1;N=500;dt=T/N;
dW=zeros(1,N);
% 布朗增量存放位置
W=zeros(1,N);
% 预分配,提高效率
dW(1)=sqrt(dt)*randn; % 循环前的初始化 W(1)=dW(1); %Matlab中数组下标从1开始,故 W(0)=0不
plot([0:dt:T],[0,W],’r-’) %绘图 xlabel(’t’,’FontSize’,16) ylabel(’W(t)’,’FontSize’,16,’Rotation’,0)
1.2 随机积分
随机积分分为Itó型随机积分和Stratonovich型随机积分。以
下假设Wiener过程 W (t),t 0 定义在概率空间 (, F , P )上,
S (t )
物理上理解,布朗运动的起因是液体的所有分子都处在运动 中,而且相互碰撞,从而微粒周围有大量的分子以微小但起伏不
定的力共同作用于它,使它被迫作不规则运动。如果用 X t 表示
微粒在时刻 t 所处位置的一个坐标,由于液体是均匀的,自然设
想 之从和时,刻因而t1 到根据t2中的心位极移限X定t2理,X可t1以是合许理多的几假乎定完X全t独2 立的X 小t1 服位从移
f (t, x) f (t, y) g(t, x) g(t, y) L2 x y ,x R, 且有E y0 2 , 则方程 (6)存在唯一解且E y(t ) 2 。
定义 2.1 (强收敛性) 若存在常数 C 0 (与 h 独立), 0 ,使得 E( y(tn ) yn ) Chp , h (0, ),
正态分布,而且对于不同时间段的位移应该是相互独立的。因此 ,布朗运动有如下定义:
定义1.1 一个随机过程 {W (t),t 0} ,它在一个微小时间间隔
t 之间内的变化为 W 。如果
1) W (0) 0;
2) W N (0, 2t ) ,其中 0为一常数。
3)对于任何两个不同时间间隔, W 的值相互独立,即独立增量。
0 t0 t1 t2 tn t ,
令tk
t k tk1(1
k
n),
max
1 k n
tk
,
若随机变量序列
n
X (tk1 )(W (tk ) W (tk1 )), n 1, 2, 3
(4)
k 1
均方收敛于唯一极限,则称
n
t
lim
n
k 1
X
(tk1 )(W
(tk
)
W
(tk1 ))
2.1 随机Taylor展开
方便起见,对如下的标量自治型随机微分方程进行讨论:
dy(t ) f ( y(t ))dt g( y(t ))dW (t ),
(10)
其中 t [t0 ,T ], X (t0 ) X 0 , X 0 R, W (t ) 是标准Wiener过
程。
随机Taylor展开式是随机微分方程数值算法的基础,Euler算
f , g 均为 [t0 ,T ]上的Borel可测函数,分别被称为漂移系数和扩散
系数。
方程(6)的积分形式为:
t
t
y(t) y(t0 )
f (s, y(s))ds
t0
g(s, y(s))dW (s),
t0
(7)
其中的随机积分为Itó型随机积分。
若将Itó型随机积分替换为Stratonovich型随机积分,则(7)式
dy(t ) f (t, y(t ))dt g(t, y(t ))dW (t ), (6)
其中
y(t0 ) y0 , t [t0 ,T ],W (t ) (W1(t ),W2(t ), ,Wd (t )) ,
E y0 2 , f : Rm [t0 ,T ] Rm , g : Rm [t0 ,T ] Rmd ,
变为
y(t) y(t0 )
t
f (s, y(s))ds
t0
t
g(s, y(s))
t0
dW (s),
(8)
对应的微分方程为
dy(t ) f (t, y(t ))dt g(t, y(t )) dW (t ), (9)
方程(9)即为Stratonovich型随机微分方程。
注:1)Itó型随机微分方程(6)和Stratonovich型随机微分方程(9)
若记随机变量 N (0,1), 则有 W t . 形式上看,当
t 0时,如同普通微积分中的情形,有:
dW dt ,
由于布朗运动是处处不可微的,此处的 dW只能视为一种简单记 法。
布朗运动的模拟
以下对一维的布朗运动进行随机模拟。一维的布朗运动可以
看 在做直质线点上在的直位线置上。作利简用单Ma随tl机ab游模动拟,布则朗W运(动t)表的示程质序点代在码时如刻下t :时
{Ft , t
0}
为 F的上升滤子(即 11
Ft
F ,且对 0
t1
t2 , Ft1
Ft2 )
,对任意 0 s t ,W ( s)关于 Ft 可测,且满足
E(W (t ) | Fs ) W (s) a.s.,
E(W (t ) W (s) | Fs ) 0 a.s., 此外,对随机过程{ X (t ), t 0},T 0, 引入以下三个条件:
X (t )dW (t )
0
(5)
为 { X (t ), t 0}关于{W (t ), t 0}在[0, t]上的Itó积分。上述定
义中,作和式(4)时不能像通常积分那样,tk 在[tk1, tk ] 中任取
,否则可能导致均方极限不存在。t (5)中取的是的[tk1, tk ]的左
端点 tk1 ,得到Itó型随机积分
h2
f
(
y(tn
))
f
(
y(tn
))
1 2
g
2
f
(
y(tn
Байду номын сангаас))
R.
(12) 求解方程(10)的Euler方法和Milstein方法均是在(12)的基础上进行 截断而得到的。
1
0.5
W(t)
0
-0.5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
t
图1 布朗运动
还可以如下进行模拟:
randn('state',100) T=1;N=500;dt=T/N;
dW=sqrt(dt)*randn(1,N); %向量化,提高运算效率
W=cumsum(dW); %累加和计算命令, W(j)=dW(1)+dW(2)+…+dW(j);j=1,…N
1.随机微分方程概述
1.1 布朗运动介绍
布朗运动是历史上最早被认真研究过的随机过程。1827 年, 英国生物学家布朗(Robert Brown)首先观察和研究了悬浮在液 体中的细小花粉微粒受到水分子连续撞击形成的运动情况,布 朗运动也因此而得名。1905 年爱因斯坦(Einstein)对它做出了合 理的物理解释并求出了微粒的转移密度,1918 年维纳(Norbert Wiener)在数学上严格地定义了布朗运动(因此它有时也称为维 纳过程)。现在布朗运动已经成为了描述随机现象的基石。