实验报告七常微分方程初值问题的数值解法
常微分方程初值问题数值解法
常微分方程初值问题的数值解法在自然科学、工程技术、经济和医学等领域中，常常会遇到一阶常微分方程初值问题：(,),,(),y f x y a x b y a y '=≤≤⎧⎨=⎩ (1) 此处f 为,x y 的已知函数，0y 是给定的初始值。
本章讨论该问题的数值解法，要求f 在区域{(,)|,}G x y a x b y =≤≤<∞内连续，并对y 满足Lipschitz 条件，从而初值问题(1)有唯一的连续可微解()y y x =，且它是适定的。
1 几个简单的数值积分法1.1 Euler 方法（1）向前Euler 公式（显式Euler 公式）10(,),0,1,2,,(),n n n n y y hf x y n y y a +=+=⎧⎨=⎩(2) 其中h 为步长。
由此便可由初值0y 逐步算出一阶常微分方程初值问题(1)的解()y y x =在节点12,,x x 处的近似值12,,y y 。
该公式的局部截断误差为2()O h ，是一阶方法。
（2）向后Euler 公式（隐式Euler 公式）1110(,),0,1,2,,(),n n n n y y hf x y n y y a +++=+=⎧⎨=⎩(3) 这是一个隐格式，也是一阶方法。
这类隐格式的计算比显格式困难，一般采用迭代法求解。
首先用向前Euler 公式提供迭代初值，然后迭代计算：(0)1(1)()111(,),(,),0,1,2,n n n n k k n n n n y y hf x y y y hf x y k +++++⎧=+⎨=+=⎩ (4)1.2 梯形方法1110[(,)(,)],2(),(0,1,2,)n n n n n n h y y f x y f x y y y a n +++⎧=++⎪⎨⎪=⎩= (5) 这也是一个隐格式，是二阶方法。
一般也采用迭代法求解。
迭代公式如下：(0)1(1)()111(,),[(,)(,)],0,1,2,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 k +++++⎧=+⎪⎨=++=⎪⎩ (6)1.3 改进的Euler 方法11110(,),[(,)(,)],0,1,2,,2(),n n n n n n n n n n y y hf x y h y y f x y f x y n y y a ++++⎧=+⎪⎪=++=⎨⎪=⎪⎩(7) 为了便于上机编程计算，(7)可改写为110(,),(,),0,1,2,,1(),2(),p n n n cn n p n p c y y hf x y y y hf x y n y y y y y a ++=+⎧⎪=+⎪⎪=⎨=+⎪⎪=⎪⎩(8) 该格式是显式，也是二阶方法。
常微分方程的数值解法实验报告
常微分方程的数值解法专业班级：信息软件 姓名：吴中原 学号：120108010002 一、实验目的1、熟悉各种初值问题的算法，编出算法程序；2、明确各种算法的精度与所选步长有密切关系；通过计算更加了解各种 算法的优越性。
二、实验题目1、根据初值问题数值算法，分别选择二个初值问题编程计算；2、试分别取不同步长，考察某节点j x处数值解的误差变化情况； 3、试用不同算法求解某初值问题，结果有何异常； 4、分析各个算法的优缺点。
三、实验原理与理论基础（一） 欧拉法算法设计对常微分方程初始问题(6-1)(6-2)用数值方法求解时，我们总是认为(6-1)、(6-2)的解存在且唯一。
欧拉法是解初值问题的最简单的数值方法。
从(6-2)式由于y (x 0) = y 0已给定，因而可以算出),()('000y x f x y =。
设x 1 = h 充分小，则近似地有：),()(')()(00001y x f x y hx y x y =≈-(6-3)记 ,n ,,i x y y i i 10 )(== 从而我们可以取),(0001y x hf y y ==作为)(1x y 的近似值。
利用1y 及f (x 1, y 1)又可以算出)(2x y 的近似值：),(1112y x hf y y +=一般地，在任意点()h n x n 11+=+处)(x y 的近似值由下式给出),(1n n n n y x hf y y +=+(6-4)这就是欧拉法的计算公式，h 称为步长。
⎪⎩⎪⎨⎧==)( ),(d d 00y x y y x f x y（二）四阶龙格-库塔法算法设计：欧拉公式可以改写为：()111,i i i i y y k k hf x y +=+⎧⎪⎨=⎪⎩，它每一步计算(),f x y 的值一次，截断误差为()2o h 。
改进的欧拉公式可以改写为：()()()11212112,,i i i i i i y y k k k hf x y k hf x h y k +⎧=++⎪⎪=⎨⎪=++⎪⎩，它每一步要计算(),f x y 的值两次，截断误差为()3o h 。
第7章 常微分方程初值问题的数值解法
例1 函数f ( t , y ) = t y 在区域D0 = {( t , y ) | 1 ≤ t ≤ 2, −3 ≤ y ≤ 4}
关于y满足Lipschitz条件，相应的Lipschitz常数可取为L = 2
3 存在性定理 定理1 设函数f ( t , y )在凸集D ⊂ R 2中有定义，若存在常数
(7.2.7)
称为显式Runge-Kutta(龙格-库塔 )方法，简称R-K方法，
其中正整数N 称为R-K方法的级，所有ci , ai , bij 都是待定 常数。
根据定义(7.2.7)，N 级R-K方法(7.9)的局部截断误差为
Rn+1 = y( t n+1 ) − y( t n ) − h∑ ci ki
dy 其斜率为 = f ( t0 , y0 ) dt ( t0 , y0 ) 由 点 斜 式 写 出 切线 方 程 dy y = y0 + ( t − t0 ) = y0 + ( t − t0 ) f ( t0 , y0 ) dt ( x0 , y )
0
等步长为h,则t1 - t0 = h, 可由切线算出 y1 ： 则 y1 = y0 + hf ( t0 , y0 ) 按此逐步计算y( tn ), 在tn +1处的值 ： yn+1 = yn + hf ( tn , yn ) y 注意： 这是“ 注意 ： 这是 “ 折 yN 线法” 而非“ 线法 ” 而非 “ 切 线法” 线法 ” 除第一个 点是曲线切线外， 点是曲线切线外 ， 其他点不是切线 y2 而是折线(如右 y1 y0 图所示)。 图所示 。
பைடு நூலகம்
则称数值解法(7.5)为显式方法。否则，称数值解法(7.3) 为隐式方法。
常微分方程初值问题数值解法
数值解法的必要性
实际应用需求
许多实际问题需要求解常微分方程初值问题，如物理、 化学、生物、工程等领域。
解析解的局限性
对于复杂问题，解析解难以求得或不存在，因此需要 采用数值方法近似求解。
数值解法的优势
未来发展的方向与挑战
高精度算法
研究和发展更高精度的算法，以提高数值解的准确性和稳定性。
并行计算
利用并行计算技术，提高计算效率，处理大规模问题。
自适应方法
研究自适应算法，根据问题特性自动调整计算精度和步长。
计算机技术的发展对数值解法的影响
1 2
硬件升级
计算机硬件的升级为数值解法提供了更强大的计 算能力。
它首先使用预估方法（如欧拉方法）得到一个 初步解，然后使用校正方法（如龙格-库塔方法） 对初步解进行修正，以提高精度。
预估校正方法的优点是精度较高，且计算量相 对较小，适用于各种复杂问题。
步长与误差控制
01
在离散化过程中，步长是一个重要的参数，它决定 了离散化的精度和计算量。
02
误差控制是数值逼近的一个重要环节，它通过设定 误差阈值来控制计算的精度和稳定性。
能够给出近似解的近似值，方便快捷，适用范围广。
数值解法的历史与发展
早期发展
早在17世纪，科学家就开始尝 试用数值方法求解常微分方程。
重要进展
随着计算机技术的发展，数值 解法在20世纪取得了重要进展， 如欧拉法、龙格-库塔法等。
当前研究热点
目前，常微分方程初值问题的 数值解法仍有许多研究热点和 挑战，如高精度算法、并行计
软件优化
软件技术的发展为数值解法提供了更多的优化手 段和工具。
微分方程初值问题的数值解法
积分法:
yk 1 yk h f ( xk , yk ) y ( x0 ) y0
积分项利用矩形公式计算
(1) y( xk 1 ) y( xk )
xk 1
xk
f (t , y(t ))dt
(★)
xk 1
xk
f (t , y(t ))dt h f ( xk , yk ) y( xk 1 ) y( xk ) h f ( xk , yk )
引言
初值问题的数值解法:求初值问题的解在一系列节点的值 y ( xn )的近似值 yn 的方法.本章数值解法的特点:都是采用“步进 式”,即求解过程顺着节点排列的次序一步步向前推进. 常微分方程初值问题: dy f ( x, y ), x [a, b] dx y ( x0 ) y0
替 f (x , y)关于 y 满足Lipschitz条件. 除了要保证(1)有唯一解外,还需保证微分方程本身是稳定的,即 (1)的解连续依赖于初始值和函数 f (x , y). 也就是说, 当初始值 y0 及函数 f (x , y)有微小变化时, 只能引起解的微小变化.
注: 如无特别说明,总假设(1)的解存在唯一且足够光滑. 在 f 连续有界, 则 f (x , y)对变量 y 可微的情形下, 若偏导数 y 可取L为
也称折线法 x
2. 梯形法
若采用梯形公式计算(★)中的积分项,则有 h y ( xk 1 ) y ( xk ) [ f ( xk , y ( xk )) f ( xk 1 , y ( xk 1 ))] 2 h yk 1 yk [ f ( xk , yk ) f ( xk 1 , yk 1 )] 2 称之为梯形公式.这是一个隐式公式,通常用迭代法求解.具体做 法: (0) (0) 先用Euler法求出初值 yk ,1 即 ,将其代入梯形公式 yk 1 yk h f ( xk , yk ) 的右端,使之转化为显式公式,即 h ( l 1) (l ) yk 1 yk [ f ( xk , yk ) f ( xk 1 , yk (☆ ) 1 )] 2
第5次实验报告(常微分方程初值问题的数值解法)
班级： 学号： 姓名： 成绩：实验5 常微分方程初值问题的数值解法实验1实验目的1）熟悉欧拉法、改进欧拉法和龙格-库塔法的原理。
2）根据以上方法，编程求解常微分方程初值问题。
2 实验内容（1）编写程序，用以上各种方法求解教材P232例7-1、习题6、11的初值问题。
(2) 使用系统自带的函数dsolve 和ode45求例7-1的符号解析解和数值解。
3实验原理求解微分方程初值问题00(,)()y f x y y x y '=⎧⎨=⎩ （1） 欧拉法（显式）：10(,)n n n n n y y hf x y x x nh +=+⎧⎨=+⎩（2） 改进欧拉法(0)1(0)111(,)(,)(,)2n n n n n n n n n n y y hf x y h y y f x y f x y ++++⎧=+⎪⎨⎡⎤=++⎪⎣⎦⎩ （3） 经典龙格-库塔法（四阶）11234121324300(22)6(,)(,)22(,)22(,)()i i i i i i i i i i 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 y y x +⎧=++++⎪⎪=⎪⎪=++⎪⎨⎪=++⎪⎪⎪=++⎪=⎪⎩4实验步骤1)建立函数文件，根据各公式编写程序；2)上机调试程序，运行程序进行计算，记录计算结果；3)分析各公式计算结果，比较各公式的优缺点。
5 程序设计欧拉法改进欧拉法function Euler1(x0,y0,h,n)%(x0,y0)：方程的初值%h：步长%n：计算的步数for i=1:nx=x0+h;yp=y0+h*f(x0,y0);yc=y0+h*f(x,yp);y=(yp+yc)/2;x %在屏幕显示每一步的x值y %在屏幕显示每一步计算的方程的数值解 x0=x;y0=y;end经典龙格-库塔法1）函数function f=f(x,y) f=y-2*x/y;end6总结注：若要更改matlab计算的数值类型，可以通过在matlab中设置实现：File -> Preferences ->Array Editor窗口中，Format 下方将Default array format设置为：long解微分方程的MATLAB命令MATLAB中主要用dsolve求微分方程的符号解析解，ode45求数值解。
常微分方程初值问题的数值解法
第七章 常微分方程初值问题的数值解法--------学习小结一、本章学习体会通过本章的学习，我了解了常微分方程初值问题的计算方法，对于解决那些很难求解出解析表达式的，甚至有解析表达式但是解不出具体的值的常微分方程非常有用。
在这一章里求解常微分方程的基本思想是将初值问题进行离散化，然后进行迭代求解。
在这里将初值问题离散化的方法有三种，分别是差商代替导数的方法、Taylor 级数法和数值积分法。
常微分方程初值问题的数值解法的分类有显示方法和隐式方法，或者可以分为单步法和多步法。
在这里单步法是指计算第n+1个y 的值时,只用到前一步的值，而多步法则是指计算第n+1个y 的值时，用到了前几步的值。
通过对本章的学习，已经能熟练掌握如何用Taylor 级数法去求解单步法中各方法的公式和截断误差，但是对线性多步法的求解理解不怎么透切，特别是计算过程较复杂的推理。
在本章的学习过程中还遇到不少问题，比如本章知识点多，公式多，在做题时容易混淆，其次对几种R-K 公式的理解不够透彻，处理一个实际问题时，不知道选取哪一种公式，通过课本里面几种方法的计算比较得知其误差并不一样，，这个还需要自己在往后的实际应用中多多实践留意并总结。
二、本章知识梳理常微分方程初值问题的数值解法一般概念步长h ，取节点0,(0,1,...,)n t t nh n M =+=，且M t T ≤，则初值问题000'(,),()y f t y t t Ty t y =≤≤⎧⎨=⎩的数值解法的一般形式是1(,,,...,,)0,(0,1,...,)n n n n k F t y y y h n M k ++==-@显示单步法7.2.1 显示单步法的一般形式1(,,),(0,1,...,1)n n n n y y h t y h n M ϕ+=+=-定理7.2.1 设增量函数(,,)n n t y h ϕ在区域00{(,,)|,||,0}D t y h t t T y h h =≤≤<∞≤≤内对变量y 满足Lipschitz 条件，即存在常数K ，使对D 内任何两点1(,,)t u h 和2(,,)t u h ，不等式1212|(,,)(,,)|||t u h t u h K u u ϕϕ-≤-成立，那么，若单步法的局部截断误差1n R +与1(1)p h p +≥同阶，即11()p n R O h ++=，则单步法的整体截断误差1n ε+与p h 同阶，即1()p n O h ε+=。
数值分析常微分方程初值问题的数值解
2、分别取步长 ，用显式欧拉法和隐式欧拉法求解初值问题
由结果分析算法的稳定性。(解析解为 )
解：当取步长 ，计算得到以下结果
准确值
显示欧拉法
隐式欧拉法
0
100.0000
100
0.5
0.0000
5776.5
0.0004
当取步长为 ，计算结果如下：
准确值
显示欧拉法
隐式欧拉法
0.0
100
100
100
0.01
60.6531
50
66.6667
0.02
36.7879
25
44.4445
0.03
22.3130
12.5
29.6296
0.04
13.5335
6.25
19.7531
0.05
8.2085
0.025
28.6505
23.7305
32.768
0.03
22.3130
17.79797
26.2144
由结果可知，当 时，显式欧拉法不太稳定，而隐式欧拉法发散；
当 时，显式欧拉法和隐式欧拉法都是稳定的。
3、选择某常微分方程初值问题的数值方法计算 的近似值，并保证有四位有效数字。
, ,
t
0.1
0.2
0.3
0.15
0.0553
-337.5
2.3324
0.2
0.0045506.20.66640.250.0004
-759.4
0.1904
0.3
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
实验报告七常微分方程初值问题的数值解法 Document number【SA80SAB-SAA9SYT-SAATC-SA6UT-SA18】浙江大学城市学院实验报告课程名称 数值计算方法实验项目名称 常微分方程初值问题的数值解法实验成绩 指导老师（签名 ） 日期 2015/12/16一. 实验目的和要求1．用Matlab 软件掌握求微分方程数值解的欧拉方法和龙格－库塔方法； 2． 通过实例学习用微分方程模型解决简化的实际问题。
二. 实验内容和原理编程题2-1要求写出Matlab 源程序(m 文件)，并有适当的注释语句；分析应用题2-2，2-3，2-4，2-5要求将问题的分析过程、Matlab 源程序和运行结果和结果的解释、算法的分析写在实验报告上。
2-1 编程编写用向前欧拉公式和改进欧拉公式求微分方程数值解的Matlab 程序，问题如下： 在区间[],a b 内(1)N +个等距点处，逼近下列初值问题的解，并对程序的每一句添上注释语句。
Euler 法 y=euler(a,b,n,y0,f,f1,b1)改进Euler 法 y=eulerpro(a,b,n,y0,f,f1,b1)2-2 分析应用题假设等分区间数100n =，用欧拉法和改进欧拉法在区间[0,10]t ∈内求解初值问题()()20(0)10y t y t y '=-⎧⎨=⎩并作出解的曲线图形，同时将方程的解析解也画在同一张图上，并作比较，分析这两种方法的精度。
2-3 分析应用题用以下三种不同的方法求下述微分方程的数值解，取10h =画出解的图形，与精确值比较并进行分析。
1）欧拉法；2）改进欧拉法；3）龙格－库塔方法；2-4 分析应用题考虑一个涉及到社会上与众不同的人的繁衍问题模型。
假设在时刻t (单位为年)，社会上有人口()x t 人，又假设所有与众不同的人与别的与众不同的人结婚后所生后代也是与众不同的人。
而固定比例为r 的所有其他的后代也是与众不同的人。
如果对所有人来说出生率假定为常数b ，又如果普通的人和与众不同的人的婚配是任意的，则此问题可以用微分方程表示为：其中变量()()()i p t x t x t =表示在时刻t 社会上与众不同的人的比例，()i x t 表示在时刻t 人口中与众不同的人的数量。
1）假定(0)0.01,0.02p b ==和0.1r =，当步长为1h =年时，求从0t =到50t =解()p t 的近似值，并作出近似解的曲线图形。
2）精确求出微分方程的解()p t ，并将你当50t =时在分题(b)中得到的结果与此时的精确值进行比较。
【MATLAB 相关函数】求微分方程的解析解及其数值的代入dsolve(‘egn1’， ‘egn2’， ‘x ’)subs (expr, {x,y,…}, {x1,y1,…} )其中‘egn i ’表示第i 个方程，‘x ’表示微分方程中的自变量，默认时自变量为t 。
subs 命令中的expr 、x 、y 为符合型表达式，x 、y 分别用数值x1、x2代入。
>> syms x y z>> subs('x+y+z',{x,y,z},{1,2,3})ans =6>> syms x>> subs('x^2',x,2)ans =4>> s=dsolve(‘12Dy y ∧=+’, ‘(0)1y =’, ‘x ’)ans =>> syms x>> subs(s,x,2)ans =右端函数(,)f x y 的自动生成f= inline(‘expr ’, ’var1’, ‘var2’,……)其中’expr ’表示函数的表达式，’var1’, ‘var2’ 表示函数表达式中的变量，运行该函数，生成一个新的函数表达式为f (var1, var2, ……)。
>> f=inline('x+3*y','x','y')f =Inline function:f(x,y) = x+3*y>> f(2,3)ans =114，5阶龙格－库塔方法求解微分方程数值解[t,x]=ode45(f,ts,x0,options)其中f 是由待解方程写成的m 文件名；x0为函数的初值；t,x 分别为输出的自变量和函数值(列向量)，t 的步长是程序根据误差限自动选定的。
若ts=[t0,t1,t2,…,tf]，则输出在自变量指定值，等步长时用ts=t0:k:tf ，输出在等分点；options 用于设定误差限(可以缺省，缺省时设定为相对误差310-，绝对误差610-)，程序为：options=odeset(‘reltol ’,rt,’abstol ’,at)，这里rt,at例：在命令窗口执行>> odefun = inline (‘2*y t y -’, ‘t ’, ‘y ’)；>> [],45(,[0,4],1)t y ode odefun =；ans =>> plot(t ,y ,‘o-’,) %解函数图形表示>> 45(,[0,4],1)ode odefun %不用输出变量，则直接输出图形>> [],45(,0:4,1)t y ode odefun =；[],t yans =三. 操作方法与实验步骤（包括实验数据记录和处理）2-1 编程编写用向前欧拉公式和改进欧拉公式求微分方程数值解的Matlab 程序，问题如下： 在区间[],a b 内(1)N +个等距点处，逼近下列初值问题的解，并对程序的每一句添上注释语句。
Euler 法 y=euler(a,b,n,y0,f,f1,b1)改进Euler 法 y=eulerpro(a,b,n,y0,f,f1,b1)Euler 法y=euler(a,b,n,y0,f,f1,b1)y=zeros(1,n+1);y(1)=y0;h=(b-a)/n;x=a:h:b;for i=1:n;y(i+1)=y(i)+h*f(x(i),y(i)); endplot(x,y)hold on% 求微分方程的精确解x1=linspace(a,b,100);'精确解为's=dsolve(f1,b1,'x')syms xy1=zeros(1,100);fori=1:100y1(i)=subs(s,x,x1(i)); endplot(x1,y1,'r')title('红色代表精确解')改进Euler法y=eulerpro(a,b,n,y0,f,f1,b1) % 求微分方程的数值解y=zeros(1,n+1);y(1)=y0;h=(b-a)/n;x=a:h:b;fori=1:n;T1=f(x(i),y(i));T2=f(x(i+1),y(i)+h*T1);y(i+1)=y(i)+(h/2)*(T1+T2); endplot(x,y)hold on% 求微分方程的精确解x1=linspace(a,b,100);'精确解为's=dsolve(f1,b1,'x')syms xy1=zeros(1,100);for i=1:100y1(i)=subs(s,x,x1(i));endplot(x1,y1,'r')title('红色代表精确解')2-2分析应用题假设等分区间数100n =，用欧拉法和改进欧拉法在区间[0,10]t ∈内求解初值问题()()20(0)10y t y t y '=-⎧⎨=⎩ 并作出解的曲线图形，同时将方程的解析解也画在同一张图上，并作比较，分析这两种方法的精度。
(1)向前欧拉法>> euler(0,10,100,10,inline('y-20','x','y'),'Dy=y-20','y(0)=10')ans =精确解为s =20 - 10*exp(x)ans =+005 *Columns 1 through 8Columns 9 through 16Columns 17 through 24Columns 25 through 32Columns 33 through 40Columns 41 through 48Columns 49 through 56Columns 57 through 64Columns 65 through 72Columns 73 through 80Columns 81 through 88Columns 89 through 96Columns 97 through 101（2）改进欧拉法>> eulerpro(0,10,100,10,inline('y-20','x','y'),'Dy=y-20','y(0)=10') ans =精确解为s =20 - 10*exp(x)ans =+005 *Columns 1 through 8Columns 9 through 16Columns 17 through 24Columns 25 through 32Columns 33 through 40Columns 41 through 48Columns 49 through 56Columns 57 through 64Columns 65 through 72Columns 73 through 80Columns 81 through 88Columns 89 through 96Columns 97 through 101改进欧拉法的精度比向前欧拉法更高。
2-3分析应用题用以下三种不同的方法求下述微分方程的数值解，取10h =画出解的图形，与精确值比较并进行分析。
1）欧拉法；2）改进欧拉法；2-4分析应用题考虑一个涉及到社会上与众不同的人的繁衍问题模型。
假设在时刻t (单位为年)，社会上有人口()x t 人，又假设所有与众不同的人与别的与众不同的人结婚后所生后代也是与众不同的人。
而固定比例为r 的所有其他的后代也是与众不同的人。
如果对所有人来说出生率假定为常数b ，又如果普通的人和与众不同的人的婚配是任意的，则此问题可以用微分方程表示为： 其中变量()()()i p t x t x t =表示在时刻t 社会上与众不同的人的比例，()i x t 表示在时刻t 人口中与众不同的人的数量。