计算方法 常微分方程初值问题数值解法Euler公式龙格库塔法
- 格式:pptx
- 大小:605.42 KB
- 文档页数:43
毕业论文文献综述信息与计算科学常微分方程初值问题的Runge-Kutta 解法一、前言部分常微分方程在很多学科领域内有着重要的作用,如自动控制、各种电子学装置的设计、弹道的计算、飞机和导弹飞行的稳定性的研究、化学反应过程稳定性的研究等等,这些问题都可以化为求微分方程的解,或者化为研究解的性质的问题。
这些问题都包含某个变量关于另一个变量的变化。
大多数这样的问题需要求解一个初值问题,即求解满足给定初值条件的微分方程[]1。
我们更多的是使用逼近原问题的解的方法来逼近原方程的解。
因逼近方法给出更精确的结果和实际的误差信息,所以更经常被使用。
一些典型的常微分初值问题的数值求解问题的方法有:单步法和线性多步法。
在求解区间[]b a ,上取定节点b x x x x x a N N =<<<<<=-1210令n n n x x x -=∆+1,称为积分网格的步长。
用N y y y ,,,10 表示初值问题精确解()x y 在节点N x x x ,,,10 上函数值()()()()N x y x y x y x y ,,,,210 的近似值。
对给定的数值积分方法,{}n y 中的各个n y 市按某一个递推算法确定的。
一个递推算法,如果在用它计算1+n y 时只用到已经求出的诸值n y y y ,,10中的n y ,而无须使用其余值110,,-n y y y 中的任何一个,则称此算法为单步方法[]2。
相反的则是多步法。
单步法主要有欧拉法和Runge-Kutta 法,多步法主要有Adams 法和Milne 法等[]3。
本文综述常微分初值问题初值问题的数值解法及其误差估计(相容性、稳定性和收敛性分析),重点介绍了Runge-Kutta 法。
二、主题部分2.1 常微分方程的初值问题概述[]114-常微分方程在微积分概念出现后即已出现。
从莱布尼兹专门研究用变量变换解决一阶微分方程的求解问题的“求通解”时代,到1841年刘维尔的里卡蒂方程不存在一般初等解和柯西的初值问题的提出,常微分方程转向“求定解”时代。
《数值分析》课程设计常微分方程初值问题数值解的实现和分析—四阶Runge-kutta方法及预估-校正算法常微分方程初值问题数值解的实现和分析—四阶Runge-kutta方法及预估-校正算法摘要求解常微分方程的初值问题,Euler方法,改进的Euler方法及梯形方法精度比较低,所以本文构造高精度单步的四级Runge-kutta方法及高精度的多步预估—校正算法及其Matlab编程来实现对常微分方程初值问题的求解,使在求解常微分方程时,对以前积分方法的收敛速度及精度都有了很高的提高。
关键词:Runge-kutta方法,Adams方法,预估—校正算法,Matlab目录1.前言 (1)2. 几个简单的数值积分法 (2)2.1Runge-kutta方法 (2)2.1.1 Runge-kutta方法的应用 (5)2.2预估—校正算法 (7)2.2.1 Adams数值积分方法简介及预估—校正算法 (7)2.2.2 预估—校正算法的应用 (12)3. 结果分析 (16)总结 (17)参考文献 (18)英文原文和中文翻译 (19)1英文原文 (19)2中文翻译 (20)1.前言常微分方程的初值问题是微分方程定解问题的一个典型代表,以下面的例子介绍常微分方程初值问题数值解的基本思想和原理。
例1.1 一重量垂直作用于弹簧所引起的震荡,当运动阻力与速度的平方成正比时,可借助如下二阶常微分方程描述若令和,则上述二阶常微分方程可化成等价的一阶常微分方程组类似于例1.1,对于m阶常微分方程其中。
若定义可得如下等价的一阶常微分方程组我们知道多数常微分方程主要靠数值解法。
所谓数值解法,就是寻求解在一系列离散节点上的近似值。
相邻两个节点之间的间距称为步长[1]。
2. 几个简单的数值积分法2.1 Runge-kutta方法Runge在1985年提出了一种基于Euler折线法的新的数值方法,此后这种新的数值方法又经过其同胞K.Heun和Kutta的努力[2],发展完善成为后世所称的Runge-kutta 方法。
《计算机数学基础》数值部分第五单元辅导14 常微分方程的数值解法一、重点内容 1. 欧拉公式:),...,,,(),()(1-210=⎩⎨⎧+=+=≈01+1+n k kh x x y x hf y y x y kk k k k k局部截断误差是O (h 2)。
2. 改进欧拉公式:预报-校正公式:⎪⎩⎪⎨⎧++=+=++++)],(),([2),(1111k k k k k k k k k k y x f y x f hy y y x hf y y 校正值预报值即 ))],(,(),([211k k k k k k k k y x hf y x f y x f hy y +++=++ 或表成平均的形式:⎪⎪⎪⎩⎪⎪⎪⎨⎧+21=+=+=1+1+)(),(),(c p k p k k c k k k p y y y y x hf y y y x hf y y改进欧拉法的局部截断误差是O (h 3)3. 龙格-库塔法二阶龙格-库塔法的局部截断误差是O (h 3) 三阶龙格-库塔法的局部截断误差是O (h 4) 四阶龙格−库塔法公式: )22(643211κκκκ++++=+hy y k k其中 κ1=f (x k ,y k );κ2=f (x n +12h ,y k +21h κ1);κ3=f (x k +12h ,y n +21h κ2);κ4=f (x k +h ,y k +h κ3)四阶龙格-库塔法的局部截断误差是O (h 5)。
二、实例例1 用欧拉法解初值问题⎩⎨⎧1=060≤≤0--='2)().(y x xy y y ,取步长h =0.2。
计算过程保留4位小数。
解h =0.2, f (x )=-y -xy 2。
首先建立欧拉迭代格式),,)((.),(210=-420=--=+=21+k y x y y hx hy y y x hf y y k k k kk k k k k k k当k =0,x 1=0.2时,已知x 0=0,y 0=1,有y (0.2)≈y 1=0.2×1(4-0×1)=0.8000当k =1,x 2=0.4时,已知x 1=0.2, y 1=0.8,有 y (0.4)≈y 2=0.2×0.8×(4-0.2×0.8)=0.614 4 当k =2,x 3=0.6时,已知x 2=0.4,y 2=0.6144,有 y (0.6)≈y 3=0.2×0.6144×(4-0.4×0.4613)=0.8000例2 用欧拉预报-校正公式求解初值问题⎩⎨⎧1=10=++'2)(sin y x y y y ,取步长h =0.2,计算y (0.2),y (0.4)的近似值,计算过程保留5位小数。
MATLAB常微分⽅程数值解——欧拉法、改进的欧拉法与四阶龙格库塔⽅法MATLAB常微分⽅程数值解作者:凯鲁嘎吉 - 博客园1.⼀阶常微分⽅程初值问题2.欧拉法3.改进的欧拉法4.四阶龙格库塔⽅法5.例题⽤欧拉法,改进的欧拉法及4阶经典Runge-Kutta⽅法在不同步长下计算初值问题。
步长分别为0.2,0.4,1.0.matlab程序:function z=f(x,y)z=-y*(1+x*y);function R_K(h)%欧拉法y=1;fprintf('欧拉法:x=%f, y=%f\n',0,1);for i=1:1/hx=(i-1)*h;K=f(x,y);y=y+h*K;fprintf('欧拉法:x=%f, y=%f\n',x+h,y);endfprintf('\n');%改进的欧拉法y=1;fprintf('改进的欧拉法:x=%f, y=%f\n',0,1);for i=1:1/hx=(i-1)*h;K1=f(x,y);K2=f(x+h,y+h*K1);y=y+(h/2)*(K1+K2);fprintf('改进的欧拉法:x=%f, y=%f\n',x+h,y);endfprintf('\n');%龙格库塔⽅法y=1;fprintf('龙格库塔法:x=%f, y=%f\n',0,1);for i=1:1/hx=(i-1)*h;K1=f(x,y);K2=f(x+h/2,y+(h/2)*K1);K3=f(x+h/2,y+(h/2)*K2);K4=f(x+h,y+h*K3);y=y+(h/6)*(K1+2*K2+2*K3+K4);fprintf('龙格库塔法:x=%f, y=%f\n',x+h,y);end结果:>> R_K(0.2)欧拉法:x=0.000000, y=1.000000欧拉法:x=0.200000, y=0.800000欧拉法:x=0.400000, y=0.614400欧拉法:x=0.600000, y=0.461321欧拉法:x=0.800000, y=0.343519欧拉法:x=1.000000, y=0.255934改进的欧拉法:x=0.000000, y=1.000000改进的欧拉法:x=0.200000, y=0.807200改进的欧拉法:x=0.400000, y=0.636118改进的欧拉法:x=0.600000, y=0.495044改进的欧拉法:x=0.800000, y=0.383419改进的欧拉法:x=1.000000, y=0.296974龙格库塔法:x=0.000000, y=1.000000龙格库塔法:x=0.200000, y=0.804636龙格库塔法:x=0.400000, y=0.631465龙格库塔法:x=0.600000, y=0.489198龙格库塔法:x=0.800000, y=0.377225龙格库塔法:x=1.000000, y=0.291009>> R_K(0.4)欧拉法:x=0.000000, y=1.000000欧拉法:x=0.400000, y=0.600000欧拉法:x=0.800000, y=0.302400改进的欧拉法:x=0.000000, y=1.000000改进的欧拉法:x=0.400000, y=0.651200改进的欧拉法:x=0.800000, y=0.405782龙格库塔法:x=0.000000, y=1.000000龙格库塔法:x=0.400000, y=0.631625龙格库塔法:x=0.800000, y=0.377556>> R_K(1)欧拉法:x=0.000000, y=1.000000欧拉法:x=1.000000, y=0.000000改进的欧拉法:x=0.000000, y=1.000000改进的欧拉法:x=1.000000, y=0.500000龙格库塔法:x=0.000000, y=1.000000龙格库塔法:x=1.000000, y=0.303395注意:在步长h为0.4时,要将for i=1:1/h改为for i=1:0.8/h。
常微分方程的数值求解在数学中,常微分方程是一类重要的数学模型,通常用来描述物理、化学、生物等自然现象中的变化规律。
对于一些复杂的微分方程,无法通过解析方法进行求解,这时候就需要借助数值方法来近似求解。
本文将介绍常微分方程的数值求解方法及其应用。
一、数值求解方法常微分方程的数值求解方法主要包括欧拉法、改进的欧拉法、龙格-库塔法等。
欧拉法是最简单也是最常用的数值求解方法,其基本思想是根据微分方程的导数近似求解下一个时间步上的解,并通过逐步迭代来得到整个解的数值近似。
改进的欧拉法在欧拉法的基础上做出了一定的修正,提高了数值求解的精度。
而龙格-库塔法则是一种更加精确的数值求解方法,通过考虑多个点的斜率来进行求解,从而减小误差。
二、应用领域常微分方程的数值求解方法在科学研究和工程实践中有着广泛的应用。
在物理学中,通过数值求解微分方程可以模拟天体运动、粒子运动等现象;在生物学领域,可以模拟生物种群的增长和变化规律;在工程领域,可以通过数值求解微分方程来设计控制系统、优化结构等。
三、实例分析以一个简单的一阶常微分方程为例:dy/dx = -y,初始条件为y(0) = 1。
我们可以用欧拉法来进行数值求解。
将时间间隔取为0.1,通过迭代计算可以得到y(1)的近似值为0.367。
而利用改进的欧拉法或者龙格-库塔法可以得到更加精确的数值近似。
这个例子展示了数值方法在解决微分方程问题上的有效性。
四、总结常微分方程是求解自然界中变化规律的重要数学工具,而数值方法则是解决一些难以解析求解的微分方程的有效途径。
通过本文的介绍,读者可以了解常微分方程的数值求解方法及其应用,希望可以对相关领域的研究和实践有所帮助。
至此,关于常微分方程的数值求解的文章正文部分结束。
浙江大学城市学院实验报告课程名称数值计算方法实验项目名称常微分方程初值问题的数值解法 实验成绩指导老师签名日期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=eulera,b,n,y0,f,f1,b1改进Euler 法y=eulerproa,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 ’ subsexpr,{x,y,…},{x1,y1,…}其中‘egn i ’表示第i 个方程,‘x ’表示微分方程中的自变量,默认时自变量为t ; subs 命令中的expr 、x 、y 为符合型表达式,x 、y 分别用数值x1、x2代入; >>symsxyz>>subs'x+y+z',{x,y,z},{1,2,3} ans= 6>>symsx>>subs'x^2',x,2 ans= 4>>s=dsolve‘12Dy y ∧=+’,‘(0)1y =’,‘x ’ ans= >>symsx >>subss,x,2 ans=右端函数(,)f x y 的自动生成f=inline ‘expr ’,’var1’,‘var2’,……其中’expr ’表示函数的表达式,’var1’,‘var2’表示函数表达式中的变量,运行该函数,生成一个新的函数表达式为fvar1,var2,……; >>f=inline'x+3y','x','y' f=Inlinefunction: fx,y=x+3y >>f2,3 ans= 114,5阶龙格-库塔方法求解微分方程数值解t,x=ode45f,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 分别为设定的相对误差和绝对误差;常用选项见下表;选项名 功能 可选值 省缺值 AbsTol 设定绝对误差正数 RelTol 设定相对误差 正数InitialStep 设定初始步长 正数 自动 MaxStep设定步长上界正数MaxOrder 设定ode15s 的最高阶数 1,2,3,4,5 5 Stats 显示计算成本统计 on,off off BDF 设定ode15s 是否用反向差分on,offoff例:在命令窗口执行>>odefun =inline ‘2*y t y -’,‘t ’,‘y ’;>>[],45(,[0,4],1)t y ode odefun =;ans=>>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=eulera,b,n,y0,f,f1,b1改进Euler 法y=eulerproa,b,n,y0,f,f1,b1Euler 法y=eulera,b,n,y0,f,f1,b1 y=zeros1,n+1; y1=y0; h=b-a/n; x=a:h:b; fori=1:n; yi+1=yi+hfxi,yi; end plotx,y holdon%求微分方程的精确解 x1=linspacea,b,100; '精确解为' s=dsolvef1,b1,'x' symsxy1=zeros1,100; for i=1:100y1i=subss,x,x1i; endplotx1,y1,'r'title'红色代表精确解'改进Euler 法y=eulerproa,b,n,y0,f,f1,b1 %求微分方程的数值解 y=zeros1,n+1; y1=y0; h=b-a/n; x=a:h:b; fori=1:n; T1=fxi,yi; T2=fxi+1,yi+hT1; yi+1=yi+h/2T1+T2; end plotx,y holdon%求微分方程的精确解 x1=linspacea,b,100; '精确解为' s=dsolvef1,b1,'x' symsxy1=zeros1,100; fori=1:100 y1i=subss,x,x1i; endplotx1,y1,'r'title'红色代表精确解' 2-2分析应用题假设等分区间数100n =,用欧拉法和改进欧拉法在区间[0,10]t ∈内求解初值问题()()20(0)10y t y t y '=-⎧⎨=⎩并作出解的曲线图形,同时将方程的解析解也画在同一张图上,并作比较,分析这两种方法的精度;1向前欧拉法>>euler0,10,100,10,inline'y-20','x','y','Dy=y-20','y0=10' ans= 精确解为 s= 20-10expx ans= +005Columns1through8(2)改进欧拉法>>eulerpro0,10,100,10,inline'y-20','x','y','Dy=y-20','y0=10' ans= 精确解为 s= 20-10expx ans= +005Columns1through8改进欧拉法的精度比向前欧拉法更高; 2-3分析应用题用以下三种不同的方法求下述微分方程的数值解,取10h = 画出解的图形,与精确值比较并进行分析; 1欧拉法; 2改进欧拉法;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 中得到的结果与此时的精确值进行比较;1>>euler0,50,50,,inline'','t','p','Dp=','p0= 1' ans= 精确解为 s=1-99/100expx/500 ans=Columns1through82>>dsolve'Dp=','p0=','t' ans=1-99/100expt/500 >>1-99/100exp ans=与欧拉法求得的精确值差0,0001四.实验结果与分析。
微分方程的常用数值解法摘要:微分方程是数学中的一种重要的方程类型,它能描述自然现象和工程问题中的许多变化规律。
但是大多数微分方程解法是无法用解析的方式求解的,因此需要借助数值解法来近似求解。
本文将介绍微分方程的常用数值解法。
关键词:欧拉方法;龙格-库塔方法;微分方程;常用数值解法一、微分方程数值解方法微分方程数值解法是数学中的重要部分。
欧拉方法、龙格-库塔方法和二阶龙格-库塔方法是常用的微分方程数值解法,下面就分别介绍这三种方法。
(一)欧拉方法欧拉方法是解初值问题的一种简单方法,它是欧拉用的第一种数值方法,也叫向前欧拉法。
欧拉方法是利用微分方程的定义式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值。
第8章 常微分方程数值解法本章主要内容:1.欧拉法、改进欧拉法. 2.龙格-库塔法。
3.单步法的收敛性与稳定性。
重点、难点一、微分方程的数值解法在工程技术或自然科学中,我们会遇到的许多微分方程的问题,而我们只能对其中具有较简单形式的微分方程才能够求出它们的精确解。
对于大量的微分方程问题我们需要考虑求它们的满足一定精度要求的近似解的方法,称为微分方程的数值解法。
本章我们主要讨论常微分方程初值问题⎪⎩⎪⎨⎧==00)(),(yx y y x f dx dy的数值解法。
数值解法的基本思想是:在常微分方程初值问题解的存在区间[a,b]内,取n+1个节点a=x 0<x 1<…<x N =b (其中差h n = x n –x n-1称为步长,一般取h 为常数,即等步长),在这些节点上把常微分方程的初值问题离散化为差分方程的相应问题,再求出这些点的上的差分方程值作为相应的微分方程的近似值(满足精度要求)。
二、欧拉法与改进欧拉法欧拉法与改进欧拉法是用数值积分方法对微分方程进行离散化的一种方法。
将常微分方程),(y x f y ='变为()*+=⎰++11))(,()()(n xn x n n dtt y t f x y x y1.欧拉法(欧拉折线法)欧拉法是求解常微分方程初值问题的一种最简单的数值解法。
欧拉法的基本思想:用左矩阵公式计算(*)式右端积分,则得欧拉法的计算公式为:Nab h N n y x hf y y n n n n -=-=+=+)1,...,1,0(),(1 欧拉法局部截断误差11121)(2++++≤≤''=n n n n n x x y h R ξξ或简记为O (h 2)。
我们在计算时应注意欧拉法是一阶方法,计算误差较大。
欧拉法的几何意义:过点A 0(x 0,y 0),A 1(x 1,y 1),…,A n (x n ,y n ),斜率分别为f (x 0,y 0),f (x 1,y 1),…,f (x n ,y n )所连接的一条折线,所以欧拉法亦称为欧拉折线法。
常微分方程初值问题的数值解法在实际应用中,对于某些微分方程,我们并不能直接给出其解析解,需要通过数值方法来求得其近似解,以便更好地理解和掌握现象的本质。
常微分方程初值问题(IVP)即为一种最常见的微分方程求解问题,其求解方法有多种,本文将对常微分方程初值问题的数值解法进行较为详细的介绍。
一、欧拉法欧拉法是最基本的一种数值解法,它采用泰勒级数展开并截断低阶项,从而获得一个差分方程近似求解。
具体来讲,设 t 为独立变量,y(t) 为函数 y 关于 t 的函数,方程为:$$y'(t) = f(t, y(t)), \qquad y(t_0) = y_0$$其中 f(t,y(t)) 为已知的函数,y(t_0) 为已知的初值。
将函数 y(t) 进行泰勒级数展开:$$y(t+h) = y(t) + hf(t, y(t)) + O(h^2)$$其中 h 表示步长,O(h^2) 表示其他高阶项。
为了使误差较小,一般取步长 h 尽可能小,于是我们可以用欧拉公式表示数值解:$$y_{n+1} = y_n + hf(t_n, y_n), \qquad y_0 = y(t_0)$$欧拉法的优点是容易理解和实现,但是由于截取低阶项且使用的单步法,所以误差较大,精度较低,在具体应用时需要慎重考虑。
二、龙格-库塔法龙格-库塔法(Runge-Kutta method)是一种多步法,比欧拉法更加精确。
龙格-库塔法的主要思想是使用不同的插值多项式来计算近似解,并且将时间步长分解,每次计算需要多次求解。
以下简要介绍二阶和四阶龙格-库塔法。
二阶龙格-库塔法将时间步长 h 分解成两步 h/2,得到近似解表达式:$$\begin{aligned} k_1 &= hf(t_n, y_n)\\ k_2 &= hf(t_n+h/2,y_n+k_1/2)\\ y_{n+1} &= y_n+k_2+O(h^3)\\ \end{aligned}$$四阶龙格-库塔法四阶龙格-库塔法是龙格-库塔法中应用最为广泛的一种方法,其需要计算的中间值较多,但是具有更高的精度。
实验八 常微分方程初值问题数值解法一、基本题科学计算中经常遇到微分方程(组)初值问题,需要利用Euler 法,改进Euler 法,Rung-Kutta 方法求其数值解,诸如以下问题:(1) ()⎪⎩⎪⎨⎧=-='004y xy y x y 20≤<x分别取h=0.1,0.2,0.4时数值解。
初值问题的精确解245x y e -=+。
(2) ()⎩⎨⎧=--='0122y y x y 01≤≤-x用r=3的Adams 显式和预 - 校式求解取步长h=0.1,用四阶标准R-K 方法求值。
(3)()()()100010321331221==-='⎪⎩⎪⎨⎧-='-='='y y y y y y y y y 10≤≤x用改进Euler 法或四阶标准R-K 方法求解取步长0.01,计算(0.05),(0.1y y y 数值解,参考结果 123(0.15)0.9880787,(0.15)0.1493359,(0.15)0.8613125y y y ≈-≈≈。
(4)利用四阶标准R- K 方法求二阶方程初值问题的数值解(I )()()⎩⎨⎧='==+'-''10,00023y y y y y 02.0,10=≤≤h x(II)()()()⎩⎨⎧='==+'--''00,10011.02y y y y y y 1.0,10=≤≤h x(III)()()⎪⎩⎪⎨⎧='=+='00,101y y e y y x 1.0,20=≤≤h x(IV)()()⎩⎨⎧='==+''00,100sin y y y y 2.0,40=≤≤h x二、应用题1. 小型火箭初始质量为900千克,其中包括600千克燃料。
火箭竖直向上发射时燃料以15千克/秒的速率燃烧掉,由此产生30000牛顿的恒定推力.当燃料用尽时引擎关闭。