解微分方程欧拉法-R-K法及其MATLAB实例

合集下载

龙哥库塔法or欧拉法求解微分方程matlab实现

龙哥库塔法or欧拉法求解微分方程matlab实现

龙哥库塔法or欧拉法求解微分⽅程matlab实现举例:分别⽤欧拉法和龙哥库塔法求解下⾯的微分⽅程我们知道的欧拉法(Euler)"思想是⽤先前的差商近似代替倒数",直⽩⼀些的编程说法即:f(i+1)=f(i)+h*f(x,y)其中h是设定的迭代步长,若精度要求不⾼,⼀般可取0.01。

在定义区间内迭代求解即可。

龙哥库塔法⼀般⽤于⾼精度的求解,即⾼阶精度的改进欧拉法,常⽤的是四阶龙哥库塔,编程语⾔如下:y(i+1)=y(i)+h*(k1+2*K2+2*k3+k4)/6;k1=f(xi,yi)k2=f(xi+h/2,yi+h*k1/2);k3=f(xi+h/2,yi+h*k2/2);k4=f(xi+h,yi+h*k3);设置终⽌条件迭代求解。

matlab实现程序如下:%% 龙哥库塔or欧拉法求解微分⽅程t=0:0.01:3; %⾃变量范围f = [];g = [];f(1) = 0.1; %f初值g(1) = 1; %g初值h=0.001;for i=1:length(t)%% 欧拉法% f(i+1) =f(i)+h*(exp(f(i)*sin(t(i)))+g(i));% g(i+1) =g(i)+h*(exp(g(i)*cos(t(i)))+f(i));%% 龙哥库塔kf1 = exp(f(i)*sin(t(i)))+g(i);%g(i)相当于常值kf2 = exp((f(i)+kf1*h/2)*sin(t(i)+h/2))+g(i);kf3 = exp((f(i)+kf2*h/2)*sin(t(i)+h/2))+g(i);kf4 = exp((f(i)+kf3*h)*sin(t(i)+h))+g(i);f(i+1) = f(i) + h*(kf1+2*kf2+2*kf3+kf4)/6;kg1 = exp(f(i)*cos(t(i)))+f(i);%f(i)相当于常值kg2 = exp((f(i)+kg1*h/2)*cos(t(i)+h/2))+f(i);kg3 = exp((f(i)+kg2*h/2)*cos(t(i)+h/2))+f(i);kg4 = exp((f(i)+kg3*h)*cos(t(i)+h))+f(i);g(i+1) = g(i) + h*(kg1+2*kg2+2*kg3+kg4)/6;endfigure(1)plot(t,f(1:length(t)),t,g(1:length(t)));legend('f函数','g函数')。

Matlab作业龙格库塔欧拉方法解二阶微分方程精选

Matlab作业龙格库塔欧拉方法解二阶微分方程精选

? 现在求解的是一个类似的问题,在这里的单摆是一种特别 的单摆,具有均匀的质量M分布在长为2的臂状摆上。
? 使用能量法(动能定理)建立方程 T ? W
1 J? 2 ? mg? h
2
? 化简得到
d 2?
dt 2
?
7.35499cos? (重力加速度取9.80665m/s2)
计算
? 边值条件y(0)=0,y'(0)=0.
接下来进行对比
计算
? 运行第三个程序:在一幅图中显示欧拉法和RK4法,随 着截断误差的积累,欧拉法产生了较大的误差
? h=0.01
h=0.0001
总结
? 通过这两种方法计算出角度峰值y=3.141593,周期是 1.777510。
? Euler方法结构简单,但是由于截断误差,使误差较大。 ? RK4是很好的方法,很稳定,由于到五阶的时候精度并
没有相应提升,所以四阶是很常用的方法。
左平成 S14060663 储建研14-2 理论力学专业
Matlab 应用
使用Euler和Rungkutta方法解 臂状摆的能量方程
1. 背景
? 定 理
M ? J?
? 化简得到
d 2?
dt 2
?
g l
sin ?
?
0
? 这化样为在,小这于样5度比的较时容候易容解易。si实n?际上简
这是一?个解二阶常微分方程的问题。
2. 问题
? 1. 使用Euler方法
? 精度随着h的减小而更高,因为向前欧拉方法的整体截 断误差与h同阶,(因为用了泰勒公式)所以欧拉方法 的稳定区域并不大。通过减小h增加了稳定性。
h=0.01
h=0.0001

matlab软件欧拉算法教程

matlab软件欧拉算法教程

y( xn1 ) y( xn ) hK
*
寻求计算平均斜率的算法
考察欧拉法,以xn的斜率值
K1 f ( xn , yn )
作为平均斜率
考察改进的欧拉法,可以将其改写为:
1 1 yi 1 yn h K1 K 2 2 2 K1 f ( x n , yn ) K 2 f ( xn h, yn hK1 )

xn
h f ( x , y( x ))dx f ( xn , y( xn )) f ( xn1 , y( xn1 )) O( h3 ) 2
h yn1 yn [ f ( xn , yn ) f ( xn1 , yn1 )] 梯形格式 2
梯形格式是显式Euler格式与隐式Euler格式的 算术平均 Euler格式是显式算法,计算量小,但精度低 梯形格式,精度较高,但是隐式算法,需要通过 迭代过程求解,计算量大
四、四阶Runge-Kutta方法
继续上述过程,可以进一步导出四阶Runge-Kutta 格式
yn1 yn h K1 2 K 2 2 K 3 K 4 , 6 K f x , y , n n 1 h K 2 f x n 1 , yn K1 , 2 2 K f x , y h K , n 1 n 2 3 2 2 K 4 f xn1 , yn hK 3 . Biblioteka 利用数值积分求积分项的方法离散
(1)左矩形法
xn1

f ( x , y( x ))dx hf ( xn , y( xn )) O( h2 )
xn
y( xn1 ) y( xn ) hf ( xn , y( xn ))

matlab 数学实验 实验报告 欧拉公式 ROSSLER微分方程

matlab 数学实验 实验报告 欧拉公式 ROSSLER微分方程

数学实验—实验报告一、实验项目: 二、实验目的和要求1、本章将对人口变化、动物种群变迁、网络系统的可靠性分析,介绍微分方程(组)的模型建立、数值解和图形解等方法,并用MATLAB 几何直观地展示各种求解方法的求解结果。

2、利用欧拉公式求解方程三、实验题目问题一:求微分方程的解析解,并画出它们的图形, y ’=y +2x , y (0)=1, 0<x <1;问题二:用向前欧拉公式和改进的欧拉公式求方程y ’=y -2x /y , y (0)=1的数值解(0≤x ≤1 , h =0.1) 要求编写程序。

问题三:Rossler 微分方程组当固定参数b=2, c=4时,试讨论随参数a 由小到大变化(如a ∈(0,0.65])而方程解的变化情况,并且画出空间曲线图形,观察空间曲线是否形成混沌状?问题四:水的流出时间一横截面积为常数A ,高为H 的水池内盛满水,由池底一横截面积为B 的小孔放水。

设水从小孔流出的速度为v=(2gh)0.5,求在任意时刻的水面高度和将水放空所需的时间。

时间t →高度h 。

问题五:考虑相互竞争模型两种相似的群体之间为了争夺有限的同一种事物来源和生存空间而进行生死存亡竞争时,往往是竞争力较弱的种群灭亡,而竞争力较强的种群达到环境容许的最大数量假设有甲、乙两个生物种群,当它们各自生存于一个自然环境中,均服从Logistics 规律。

三、实验过程问题一:用matlab 编写代码: x=[0,1]y=dsolve('Dy=y+2*x')y=dsolve('Dy=y+2*x', 'y(0)=1', 'x')ezplot(x,y)输出:y =-2*x+exp(t)*C1 (通解)y =-2*x-2+3*exp(x) 画图:x=0:0.01:1;y =-2*x-2+3*exp(x);plot(x,y)'''()x y z y x ayz b z x c =--⎧⎪=+⎨⎪=+-⎩0.10.20.30.40.50.60.70.80.91问题二: 1、分析:解:(1)解析解法得到其精确解:(2) 向前欧拉法:1(2/)n n n n n y y h y x y +=+- (1)2/n n nh y h x y =+- 迭代公式为 n+1y 1.10.2/n n n y x y =-,其中0y =y(0)=1(3)改进欧拉法:n+1nn n n n+1n +1n n n n n n n n n n n nn2n nnnn y =y +(h /2)*[(y -2x /y )+(y -2x /y )]=y +(h /2)*[(y -2x/y )+(y +h -2(x+h )/(yh ))] =y+(h /2)*[2y2x /y2(x +h )/(y+h )]=(1+h )y /2x/y(x +h )/(y+h )h h h h +--+-- 迭代公式为 n+1y 1.10.1/0.1()/()0.005n n n n n y x y x h y h =--+++,其中0y =y(0)=12、Matlab 编码x1(1)=0;y1(1)=1;y2(1)=1;h=0.1; for k=1:10 x1(k+1)=x1(k)+h;y1(k+1)=(1-h)*y1(k)+2*h*x1(k)/y1(k);y2(k+1)=(1+h)*y2(k)+(h*h)/2-h*x1(k)/y2(k)-h*(x1(k)+h)/(y2(k)+h); end x=0:0.1:1; y=(2*x+1).^(1/2);x1=x1(1:11),y=y(1:11),y1=y1(1:11),y2=y2(1:11), plot(x,y,x1,y1,'k:',x1,y2,'r--')显示图像及结果:x1 = 0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000y = 1.0000 1.0954 1.1832 1.2649 1.3416 1.4142 1.4832 1.5492 1.6125 1.6733 1.7321y1 = 1.0000 0.9000 0.8322 0.7971 0.7926 0.8143 0.8557 0.9103 0.9731 1.0402 1.1092y2 = 1.0000 1.0959 1.1847 1.2679 1.3468 1.4222 1.4948 1.5653 1.6340 1.7016 1.768300.10.20.30.40.50.60.70.80.910.811.21.41.61.82图中,蓝色曲线是精确解,红色曲线是向前欧拉法曲线,黑色曲线是改进后欧拉法曲线问题三1、matlab 编程function r=rossler(t,x) global a; global b; global c;r=[-x(2)-x(3);x(1)+a*x(2);b+x(3)*(x(1)-c)];global a; global b; global c; b=2; c=4; t0=[0,200]; for a=0:0.02:0.65[t,x]=ode45('rossler',t0,[0,0,0]); subplot(1,2,1);plot(t,x(:,1),'r',t,x(:,2),'g',t,x(:,3),'b');title('x(ºìÉ«),y(ÂÌÉ«),z(ÀºÉ«)Ëæt±ä»¯Çé¿ö');xlabel('t'); subplot(1,2,2);plot3(x(:,1),x(:,2),x(:,3))title('Ïàͼ');xlabel('x');ylabel('y');zlabel('z'); pause end当a=0时,图像如下50100150200-1-0.8-0.6-0.4-0.200.20.40.6x(红色),y(绿色),z(篮色)随t 变化情况t-0.5相图z当a=0时,(x,y,z)收敛于(0,0.5,0.5)当a=0.12时,图像如下50100150200-1.2-1-0.8-0.6-0.4-0.200.20.40.6x(红色),y(绿色),z(篮色)随t 变化情况t-0.5相图z当a=0.28时,(x,y ,z)仍然收敛,但是收敛速度大大降低。

解微分方程欧拉法,R-K法及其MATLAB实例

解微分方程欧拉法,R-K法及其MATLAB实例

解微分方程的欧拉法,龙格-库塔法及其MATLAB简单实例欧拉方法(Euler method)用以对给定初值的常微分方程(即初值问题)求解分为前进EULER法、后退EULER法、改进的EULER法。

缺点:欧拉法简单地取切线的端点作为下一步的起点进行计算,当步数增多时,误差会因积累而越来越大。

因此欧拉格式一般不用于实际计算。

改进欧拉格式:为提高精度,需要在欧拉格式的基础上进行改进。

采用区间两端的斜率的平均值作为直线方程的斜率。

改进欧拉法的精度为二阶。

算法为:微分方程的本质特征是方程中含有导数项,数值解法的第一步就是设法消除其导数值。

对于常微分方程:x∈[a,b]y(a) = y0可以将区间[a,b]分成n段,那么方程在第xi点有y'(xi) = f(xi,y(xi)),再用向前差商近似代替导数则为:在这里,h是步长,即相邻两个结点间的距离。

因此可以根据xi点和yi点的数值计算出yi+1来:i=0,1,2,L这就是向前欧拉格式。

改进的欧拉公式:将向前欧拉公式中的导数f(xi,yi)改为微元两端导数的平均,即上式便是梯形的欧拉公式。

可见,上式是隐式格式,需要迭代求解。

为了便于求解,使用改进的欧拉公式:数值分析中,龙格-库塔法(Runge-Kutta)是用于模拟常微分方程的解的重要的一类隐式或显式迭代法。

实际上,龙格-库塔法是欧拉方法的一种推广,向前欧拉公式将导数项简单取为f(xn,yn),而改进的欧拉公式将导数项取为两端导数的平均。

龙格-库塔方法的基本思想:在区间[xn,xn+1]内多取几个点,将他们的斜率加权平均,作为导数的近似。

龙格库塔法的家族中的一个成员如此常用,以至于经常被称为“RK4”或者就是“龙格库塔法”。

令初值问题表述如下。

则,对于该问题的RK4由如下方程给出:其中这样,下一个值(y n+1)由现在的值(y n)加上时间间隔(h)和一个估算的斜率的乘积决定。

该斜率是以下斜率的加权平均:k1是时间段开始时的斜率;k2是时间段中点的斜率,通过欧拉法采用斜率k1来决定y在点tn + h/2的值;k3也是中点的斜率,但是这次采用斜率k2决定y值;k4是时间段终点的斜率,其y值用k3决定。

matlab_常微分方程数值解法

matlab_常微分方程数值解法
d2x 2x2 0
dt 2
简朴问题可以求得解析解,多数实际问题靠数值求解 。
第4页
一阶常微分方程(ODE )初值问题 : ODE :Ordinary Differential Equation
dy
f
(x,
y)
dx
x0 x xn
y(x0 ) y0
数值解法就是求y(x)在某些分立旳节点 xn 上旳近似值 yn,用以近似y(xn)
x0
y0
x1 f y(x), x dx
x0
x2 f y(x), x dx
x1
y(x1) f y(x1), x1 h
第17页
同样,在[x0,xn+1] ,积分采用矩形近似,得:
y(xn1) y0
f xn1
x0
y(x), x dx
y(xn ) f y(xn ), xn h
yn y(xn )
第5页
2、欧拉近似办法
2.1 简朴欧拉(L.Euler, 1707-1783)办法。
dy
dx
f
(y, x)
y(x0 ) y0
欧拉数值算法就是由初值通过递推求解,递推求解
就是从初值开始,后一种函数值由前一种函数值得到。核 心是构造递推公式。
y0 y1 y2 yn
第6页
i 1,2,...
第36页
没有一种算法可以有效地解决所有旳 ODE 问题,因此 MATLAB 提供了多种ODE函数。
函数 ODE类
特点
阐明

ode45
非刚性 单步法;4,5 阶 R-K 措施;合计 大部分场合旳首选措施
截断误差为 (△x)3
ode23
非刚性 单步法;2,3 阶 R-K 措施;合计 使用于精度较低旳情形

微分方程数值解之欧拉法在MATLAB下的应用

!科技风 "#"$ 年 % 月
科教论坛 !"#!$%&$'(') *+&,-./&$01$21(3$&)%)$$%%(3
微分方程数值解之欧拉法在 J7<97F下的应用
梁春叶4王桥明4孙远通4叶晓艳4曾宝莹
玉林师范学院数学与统计学院!广西玉林!(#+$$$
摘4要微分方程在实际应用中十分广泛涉及领域众多但对于微分方程的数值解的计算仍然有很大挑战 本文着重 对微分方程数值解求解的常用的一类基础方法欧拉法进行在 EH#[HC的应用下的一个简单介绍
求出$ 以下为在 50;"0)中实现!
F8/F68& 6'
+'
Copyright©博看网 . All Rights Reserved.
科教论坛
%I&Fh(J>>& 8'hh63t3i8326% KG&( h+ 8&$'hh( 27, % 8h(FG*H%& %I&F#KG&(' 由此得该方程的精确解为!
科技风 "#"$ 年 % 月
8hFJ/E*J>8& 8' 由此得解析解为!8h67 i63 N6N$' $
)% $) $) )% 使用欧拉法求解#如下! A h%3$% >hu& 6#8' 63t(i63t)% 6h+$!A!$%, % -hFJ9%& 6#)' i$% 8$ h+) 27#9%$GF&$#-' , % >G$& h$!8$& &N$'h8$& &' NA!>& 6& &' #8$& &' ' % %&( 8h63t732)%i63t332$)N632$)N$' 2)%% E*G7& 6#8#pGBp#6#8$#pN@p' *%L%&(& p,6?K7p#p,+*%$p' 运行程序可输出结果为!

matlab求解微分方程数值解与解析解

matlab求解微分方程数值解与解析解微分方程是数学中的重要内容,它描述了物理、工程、经济等领域中的许多现象和问题。

在实际应用中,我们经常需要求解微分方程的解析解或数值解。

本文将以Matlab为工具,探讨如何求解微分方程并对比解析解与数值解的差异。

一、引言微分方程是描述自然界中许多现象和问题的数学语言,它包含了未知函数及其导数与自变量之间的关系。

微分方程的求解可以帮助我们了解问题的性质和变化规律,并为实际应用提供参考。

在许多情况下,微分方程的解析解很难求得,这时我们可以利用计算机进行数值求解。

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

它通过将微分方程转化为差分方程,然后利用离散的点逼近连续的解。

具体步骤如下:(1)将微分方程转化为差分方程,即用近似的导数代替真实的导数;(2)选择初始条件,即确定初始点的值;(3)选择步长和求解区间,即确定求解的范围和步长;(4)使用迭代公式计算下一个点的值;(5)重复步骤(4),直到达到指定的求解区间。

2.改进的欧拉法欧拉法存在精度较低的问题,为了提高精度,可以使用改进的欧拉法。

改进的欧拉法是通过使用两次导数的平均值来计算下一个点的值,从而提高了数值解的精度。

3.龙格-库塔法龙格-库塔法是一种常用的数值求解微分方程的方法,它通过使用多个点的导数来逼近连续解。

龙格-库塔法的步骤如下:(1)选择初始条件和步长;(2)使用迭代公式计算下一个点的值;(3)计算下一个点的导数;(4)根据导数的值和步长计算下一个点的值;(5)重复步骤(3)和(4),直到达到指定的求解区间。

龙格-库塔法的精度较高,适用于求解一阶和高阶微分方程。

三、微分方程的解析解解析解是指能够用公式或函数表示的方程的解。

有些微分方程具有解析解,可以通过数学方法求得。

例如,一阶线性常微分方程和某些特殊类型的二阶微分方程等。

解析解的优势在于精确性和直观性,能够帮助我们深入理解问题的本质。

7-17-2微分方程01欧拉法02R-K法

yn +1 = yn + h f ( xn , yn ) = y ( xn ) + h y′( xn )
1 y ′′ (ξ ) h 2 2
二阶泰勒公式
y ( xn +1 ) = y ( xn ) + y ′ ( xn ) h +
两式相减, 两式相减,由设 yn=y(xn ) ,有
h2 y ( xn +1 ) − yn +1 = y′′ (ξ ) = O ( h 2 ) 2
一阶常微分方程初 值问题
dy = f ( x, y ) dx y ( x0 ) = y0
的数值解法。 基 本思想 :常 微分方 程初值 问题的 数值解 是求 微分方程 的 解 y ( x) 即 微 分 方 程 初 值 问 题 的 积 分 曲 线 ) 在 区 间 [ a , b] 中 ( , 点( 给 定 一 系 列 点 ( 节 点 ) x n = x n − 1 + hn ( n = 1, 2, L ) 上 的 近 似 值 y n 。 这 里 hn 为 x n −1 到 x n 的 步 长 , 且 h n > 0 。
等距, 设 x0 , x1 , x2 ,⋅ ⋅ ⋅, xn 等距,步长
y′( x ) ≈
y ( x + h ) ≈ y ( x ) + h y ′( x ) = y ( x ) + h f ( x , y )
令x=xn , x+h=xn+1 , y(xn)≈yn ,y(xn+1 ) ≈yn+1 ,
初值问题离散化为
h2 y ( xn +1 ) − yn +1 ≈ − y ''( xn ) 2

matlab实例讲解欧拉法求解微分方程

欧拉法是数值分析中常用的一种方法,用于求解常微分方程的数值解。

在MATLAB中,可以通过编写相应的代码来实现欧拉法求解微分方程。

下面我们将通过具体的实例来讲解MATLAB中如何使用欧拉法求解微分方程。

我们要了解欧拉法的基本原理。

欧拉法是一种通过迭代逼近微分方程解的方法,它基于微分方程的定义,通过离散化的方法逼近微分方程的解。

其基本思想是利用微分方程的导数定义,将微分方程以差分形式进行逼近。

具体而言,欧拉法通过将微分方程转化为差分方程的形式,然后通过迭代逼近得到微分方程的数值解。

接下来,我们通过一个具体的实例来讲解MATLAB中如何使用欧拉法求解微分方程。

假设我们要求解以下的一阶常微分方程:(1) dy/dx = x + y(2) y(0) = 1现在我们来编写MATLAB代码来实现欧拉法求解这个微分方程。

我们需要确定微分方程的迭代步长和迭代范围。

假设我们将x的范围取为0到10,步长为0.1。

接下来,我们可以编写MATLAB代码如下:```matlab欧拉法求解微分方程 dy/dx = x + y定义迭代步长和范围h = 0.1;x = 0:h:10;初始化y值y = zeros(1,length(x));y(1) = 1;使用欧拉法迭代求解for i = 1:(length(x)-1)y(i+1) = y(i) + h * (x(i) + y(i));end绘制图像plot(x,y,'-o');xlabel('x');ylabel('y');title('欧拉法求解微分方程 dy/dx = x + y');```在这段MATLAB代码中,我们首先定义了迭代的步长和范围,并初始化了微分方程的初始值y(0) = 1。

然后通过for循环使用欧拉法进行迭代求解微分方程,最后绘制出了微分方程的数值解的图像。

通过以上的实例讲解,我们可以看到,在MATLAB中使用欧拉法求解微分方程是非常简单而直观的。

  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

欧拉方法(Euler method)用以对给定初值的常微分方程(即初值问题)求解
分为前进EULER法、后退EULER法、改进的EULER法。

缺点:
欧拉法简单地取切线的端点作为下一步的起点进行计算,当步数增多时,误差会因积累而越来越大。

因此欧拉格式一般不用于实际计算。

改进欧拉格式:
为提高精度,需要在欧拉格式的基础上进行改进。

采用区间两端的斜率的平均值作为直线方程的斜率。

改进欧拉法的精度为二阶。

算法为:
微分方程的本质特征是方程中含有导数项,数值解法的第一步就是设法消除其导数值。

对于常微分方程:
x∈[a,b]
y(a) = y0
可以将区间[a,b]分成n段,那么方程在第xi点有y'(xi) = f(xi,y(xi)),再用向前差商近似代替导数则为:
在这里,h是步长,即相邻两个结点间的距离。

因此可以根据xi点和yi点的数值计算出yi+1来:
i=0,1,2,L
这就是向前欧拉格式。

改进的欧拉公式:
将向前欧拉公式中的导数f(xi,yi)改为微元两端导数的平均,即
上式便是梯形的欧拉公式。

可见,上式是隐式格式,需要迭代求解。

为了便于求解,使用改进的欧拉公式:
数值分析中,龙格-库塔法(Runge-Kutta)是用于模拟常微分方程的解的重要的一类隐式或显式迭代法。

实际上,龙格-库塔法是欧拉方法的一种推广,向前欧拉公式将导数项简单取为f(xn,yn),而改进的欧拉公式将导数项取为两端导数的平均。

龙格-库塔方法的基本思想:
在区间[xn,xn+1]内多取几个点,将他们的斜率加权平均,作为导数的近似。

龙格库塔法的家族中的一个成员如此常用,以至于经常被称为“RK4”或者就是“龙格库塔法”。

令初值问题表述如下。

则,对于该问题的RK4由如下方程给出:
其中
这样,下一个值(y n+1)由现在的值(y n)加上时间间隔(h)和一个估算的斜率的乘积决定。

该斜率是以下斜率的加权平均:
k1是时间段开始时的斜率;
k2是时间段中点的斜率,通过欧拉法采用斜率k1来决定y在点tn + h/2的值;
k3也是中点的斜率,但是这次采用斜率k2决定y值;
k4是时间段终点的斜率,其y值用k3决定。

当四个斜率取平均时,中点的斜率有更大的权值:
RK4法是四阶方法,也就是说每步的误差是h5阶,而总积累误差为h4阶。

注意上述公式对于标量或者向量函数(y可以是向量)都适用。

例子:
下面给出了数值求解该微分方程的简单程序。

其中y1,y2,y3,y4分别为向前欧拉公式,改进的欧拉公式,4级4阶龙格-库塔公式及精确解。

h=;
x=0:h:1;
y1=zeros(size(x));
y1(1)=1;
y2=zeros(size(x));
y2(1)=1;
y3=zeros(size(x));
y3(1)=1;
for i1=2:length(x)
y1(i1)=y1(i1-1)+h*(y1(i1-1)-2*x(i1-1)/y1(i1-1));
k1=y2(i1-1)-2*x(i1-1)/y2(i1-1);
k2=y2(i1-1)+h*k1-2*x(i1)/(y2(i1-1)+h*k1);
y2(i1)=y2(i1-1)+h*(k1+k2)/2;
k1=y2(i1-1)-2*x(i1-1)/y2(i1-1);
k2=y2(i1-1)+h*k1/2-2*(x(i1-1)+h/2)/(y2(i1-1)+h*k1/2);
k3=y2(i1-1)+h*k2/2-2*(x(i1-1)+h/2)/(y2(i1-1)+h*k2/2);
k4=y2(i1-1)+h*k3-2*(x(i1-1)+h)/(y2(i1-1)+h*k3);
y3(i1)=y3(i1-1)+(k1+2*k2+2*k3+k4)*h/6;
end
y4=sqrt(1+2*x);
%plot(x,y1,x,y2,x,y3,x,y4)
%legend('y1','y2','y3','y4')
plot(x,y4-y1,x,y4-y2,x,y4-y3)
legend('y1','y2','y3')。

相关文档
最新文档