5.2龙格库塔法

第五章 常微分方程的数值解法

5.2 龙格-库塔法

一、教学目标及基本要求

通过对本节课的学习，使学生掌握常微分方程、常微分方程方程组的数值解法。

二、教学内容及学时分配

本节课主要介绍常微分方程的数值解法。具体内容如下： 讲授内容：龙格库塔方法、收敛性与稳定性

三、教学重点难点

1．教学重点：龙格库塔方法、收敛性与稳定性。 2. 教学难点：收敛性与稳定性。

四、教学中应注意的问题

多媒体课堂教学为主。适当提问，加深学生对概念的理解。

五、正文

龙格－库塔方法

引言 龙格－库塔方法的基本思想

'11()()

()()()(,())n n n n y x y x y y x y x hf y h

ξξξ++-=?=+

令*(,())K f y ξξ=称为区间1[,]n n x x +上的平均斜率，只要对*K 提供一种算法，即可推导出一种计算公式。

欧拉公式只是取点n x 的斜率作为区间1[,]n n x x +的平均斜率*K ，精度自然很低。

考察改进的欧拉公式

11211

()22

n n y y h k k +=++1(,)n n k f x y =21(,)n n k f x h y hk =++

它利用n x ，1n x +两点的斜率取算术平均，1n x +处斜率通过已知信息n y 用欧拉公式预报得到。

可以考虑设法在1[,]n n x x +多取几个点的斜率值，将它们加权平均作为区间

1[,]n n x x +的平均斜率*K 。这就是龙格－库塔方法的基本思想。 1、二阶龙格－库塔方法

考察1[,]n n x x +内一点,01n p n x x hp p +=+<≤，希望通过,n n p x x +两个点的斜率

12,K K 加权平均得到*K ，即令112((1))n n y y h K K λλ+=+-+ 取1(,)n n K f x y =，如何预报n p x +处斜率2K ？

仿照改进的欧拉公式，先用欧拉公式预测()n p y x +的值n p y +：

1n p n y y phK +=+

然后用n p y +计算2(,)n p n p K f x y ++=，从而得

112121[(1)](,)(,)

n n n n n p n y y h K K K f x y K f x y phK λλ++=+-+==+

适当选取,p λ，使上述公式具有较高得精度。假定()n n y y x =，分别将12

,K K 泰勒展开：

'1(,)()n n n K f x y y x ==

212

'

''

2

(,)

(,)[(,)(,)(,)]()()()()

n p n n n x n n n n y n n n n K f x y pkK f x y ph f x y f x y f x y O h y x phy x O h +=+=+++=++

代入得

'2''31()()()()n n n n y y x hy x ph y x O h λ+=+++

按泰勒展开2'

''

31()()()()2

n n n n h y y x hy x y x O h +=+++

比较得，只要1

2

p λ=

，公式截断误差为3()O h 特别，当1,1/2p λ==，就是改进的欧拉公式， 改取1,1/2p λ==，

12n n y y hk +=+，1(,)n n k f x y =，21(,)22

n n h h

k f x y k =++

称为变形的欧拉公式，亦称中点公式。 2、三阶龙格－库塔方法

为提高精度，设除n p x +外再考察一点,01n q n x x hq p q +=+<≤≤，用三个点

,,n n p n q x x x ++的斜率123,,K K K 加权平均得出*123(1)K K K K λμλμ=--++

1123((1))n n y y h K K K λμλμ+=+--++

为预测n q x +的斜率值3K ，将12,K K 加权平均得出[,]n n q x x +上的平均斜率，从而得

12((1))n q n y y qh K K αα+=+-+

然后用n q y +计算3(,)n q n q K f x y ++= 设计的计算公式具有如下形式

1123121312[(1)](,)(,)

(,((1)))

n n n n n p n n q n y y h K K K K f x y K f x y phK K f x y qh K K λμλμαα+++=+--++==+=+-+

运用泰勒展开选择参数,,,,p q λμα 下列库塔公式是其中一种

112312112

312[4]

6

(,)(,)

2(,(2))

n n n n n n n q n h

y y K K K K f x y h

K f x y K K f x y h K K +++=+++==+=+-+

3、四阶龙格－库塔方法

继续上述过程，得出四阶方法，下列是经典格式

11234121123122

413[22]

6

(,)(,)2(,)2

(,)

n n n n n n n n n n h

y y K K K K K f x y h K f x y K h K f x

y K K f x y hK ++

+

+=++++==+=+=+

流程图P105 例P105

1211121312222413332(,)2()

2(,)()22()

22()

2(,)()22()

22()

(,)()()

n n n n n

n n n n n n n n n n n n n n n x K f x y y y h x h h K f x y K y K h y K h x h h K f x y K y K h y K x h K f x y hK y hK y hK +++==-

+=+=+-

++=+=+-

++=+=+-

+ 4、变步长龙格－库塔方法

步长越小，截断误差越小，但步长太小，计算量大，舍入误差也大，需要数值解法选择步长。步长选择考虑的问题： 1）如何衡量和检验计算结果的精度； 2）如何依据判定的精度来处理步长。

以经典龙格－库塔方法为例，先以步长h 求得一个近似值记为1h

n y +，经典龙格－库塔方法的截断误差为5()O h ，故有511()()h n n y x y CO h ++-≈

步长折半/2

511

()2(())2

h n n h y x y CO ++-≈ /2

/2/2/211

11111111()11()()()()1615

h h h h h h n n n n n n n n h

n n y x y y x y y y y y y x y δ++++++++++-≈?-≈-?=-- 区分两种情况处理步长：

1）δε>，继续将步长折半，直至δε<为止，取步长折半后的新值作为结果；

2）δε<，步长加倍直至δε>，取步长加倍前的老值作结果。

5、介绍： 亚当姆斯方法

单步法的特点是计算1n y +时，只用到n y 的值，此时，011,,......,,n n y y y y -的值均已算出，如果在计算1n y +时，除用n y 的值外，同时还用到11,...,n n k y y --+的值，则称此方法为多步法。

k 步亚当姆斯公式：

1

11k n n i n i i y y h f β-+-=-=+∑

当式中的10β-=时，称为亚当姆斯显式公式。当10β-≠时，称为亚当姆斯隐式公式。

小结

本节课主要介绍了解常微分方程的龙格库塔法。要求大家掌握低阶龙格库塔法。

作业：课后习题4-8.

matlab编的4阶龙格库塔法解微分方程的程序

matlab编的4阶龙格库塔法解微分方程的程序 2010-03-10 20:16 function varargout=saxplaxliu(varargin) clc,clear x0=0;xn=1.2;y0=1;h=0.1; [y,x]=lgkt4j(x0,xn,y0,h); n=length(x); fprintf(' i x(i) y(i)\n'); for i=1:n fprintf('%2d %4.4f %4.4f\n',i,x(i),y(i)); end function z=f(x,y) z=-2*x*y^2; function [y,x]=lgkt4j(x0,xn,y0,h) x=x0:h:xn; n=length(x); y1=x; y1(1)=y0; for i=1:n-1 K1=f(x(i),y1(i)); K2=f(x(i)+h/2,y1(i)+h/2*K1); K3=f(x(i)+h/2,y1(i)+h/2*K2); K4=f(x(i)+h,y1(i)+h*K3); y1(i+1)=y1(i)+h/6*(K1+2*K2+2*K3+K4); end y=y1; 结果： i x(i) y(i) 1 0.0000 1.0000 2 0.1000 0.9901 3 0.2000 0.9615 4 0.3000 0.9174 5 0.4000 0.8621 6 0.5000 0.8000 7 0.6000 0.7353 8 0.7000 0.6711 9 0.8000 0.6098 10 0.9000 0.5525 11 1.0000 0.5000 12 1.1000 0.4525 13 1.2000 0.4098

matlab 四阶龙格-库塔法求微分方程

Matlab 实现四阶龙格-库塔发求解微分方程 从理论上讲，只要函数在某区间上充分光滑，那么它可以展开为泰勒级数，因此在该区间上的函数值可用各阶导数值近似地表示出来，反之其各阶导数值也可用某些函数值的线性组合近似地表示出来。龙格-库塔法就是将待求函数)(t y 展开为泰勒级数，并用方程函数),(y f t 近似其各阶导数，从而迭代得到)(t y 的数值解。具体来说，四阶龙格-库塔迭代公式为 )22(6 143211k k k k h n n ++++=+y y ),(1n n t k y f = )2/,2/(12hk h t k n n ++=y f )2/,2/(23hk h t k n n ++=y f ),(33hk h t k n n ++=y f 实验内容： 已知二阶系统21x x = ，u x x x 5.02.04.0212+--= ，0)0()0(21==x x ，u 为单位阶跃信号。用四阶龙格-库塔法求数值解。分析步长对结果的影响。 实验总结： 实验报告要求简要的说明实验原理；简明扼要地总结实验内容；编制m 文件，并给出运行结果。报告格式请按实验报告模板编写。 进入matlab ， Step1：choose way1 or way2 way1）： 可以选择直接加载M 文件（函数M 文件）。 way2）： 点击new ——function ，先将shier （函数1文本文件）复制运行; 点击new ——function ，再将RK （函数2文本文件）运行; 点击new ——function ，再将finiRK （函数3文本文件）运行;

龙格库塔积分算法

龙格库塔法 龙格库塔法是常用于模拟常微分方程的解的重要的一类隐式或显式迭代法。这些技术由数学家C. Runge和M.W. Kutta于1900年左右发明。由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。 龙格库塔法是一种在工程上应用广泛的高精度单步算法，可以应用在物理、工程、控制、动力学中，如模糊控制、弹道分析以及分析光纤特性等，在系统仿真中得到广泛应用。 龙格库塔法源自于相应的泰勒级数方法，在每一插值节点用泰勒级数展开,其截断误差阶数也是，根据可省略更高阶的导数计算, 这种方法可构造任意阶数的龙格库塔法。其中4 阶龙格库塔法是最常用的一种方法。因为它相当精确、稳定、容易编程。在计算中一般不必使用高阶方法, 因为附加的计算误差可由增加精度来弥补。如果需要较高的精度, 可采取减小步长的方法即可。4 阶龙格库塔法的精度类似4 阶泰勒级数法的精度。 1、初值问题 对于一阶常微分方程的初值问题 根据常微分方程的理论可知，此初值问题的解在区间[a，b]上存在，且唯一。 2、离散化

取步长h=(b-a)/n，将区间[a , b]分成n个子区间： a=<=b 在其中任意两点的曲线段上，根据积分中值定理，一段光滑曲 线上至少有一点，它的斜率与整段曲线的平均斜率相同， 得=y’() (0<<1) 其中，= 可以将上式改写成y()=y()+h*K (2.1) 其中K为平均斜率，K=f() 公式（2.1）表明，如果能够确定平均斜率K，就可以根据（2.1）式得到y()的值。 欧拉法和龙格库塔法就是用不同方法确定不同精度的平均斜率K，从而求得y()的近似值。 3、Euler法 欧拉法虽然精度低，但它是最简单的一种显式单步法，也是龙 格库塔法的基础。 首先，令、为y() 及y()的近似值，并且令平均斜 率K=f()，即以点的斜率作为平均斜率K，便得到欧拉公式=+h* f() （3.1） 4、改进的欧拉法 此种方法是取、两点的斜率的平均值作为平均斜率K， 即K= ,其中、均为y()以及y()的近似值，就得到 改进后的欧拉公式（4.1）

控制系统数字仿真 四阶龙格库塔法

控制系统数字仿真 1.实验目的 1.掌握利用四阶龙格-库塔（Runge-Kutta）法进行控制系统数字仿真的方 法。 2.学习分析高阶系统动态性能的方法。 3.学习系统参数改变对系统性能的影响。 二、实验内容 已知系统结构如下图 若输入为单位阶跃函数，计算当超调量分别为5%，25%，和50%时K的取值（用主导极点方法估算），并根据确定的K值在计算机上进行数字仿真。 三、实验过程 1.计算K值 二阶系统单位阶跃响应的超调量 %100% =? 1.当σ%=5%时

解得 ζ=0.690 设主导极点 =ζa + a=0.69a+j0.72a 代入D （s ）= 32 1025s s s K +++=0中， 32(0.690.72)10(0.690.72)25(0.690.72)0 a j a a j a a j a K ++++++=解得K=31.3，a=-2.10 即1,2 1.45 1.52s j =-± 2. 当σ%=25%时 解得 ζ=0.403 设主导极点 =ζa + a=0.403a+j0.915a 代入D （s ）= 321025s s s K +++=0中， 32(0.4030.915)10(0.4030.915)25(0.4030.915)0 a j a a j a a j a K ++++++=解得K=59.5，a=-2.75 即1,2 1.11 2.53s j =-± 3. 当σ%=50%时 解得 ζ=0.215 设主导极点 =ζa + a=0.215a+j0.977a 代入D （s ）= 321025s s s K +++=0中， 32(0.2150.977)10(0.2150.977)25(0.2150.977)0 a j a a j a a j a K ++++++=解得K=103，a=-3.48

5.2龙格库塔法

第五章 常微分方程的数值解法 5.2 龙格-库塔法 一、教学目标及基本要求 通过对本节课的学习，使学生掌握常微分方程、常微分方程方程组的数值解法。 二、教学内容及学时分配 本节课主要介绍常微分方程的数值解法。具体内容如下： 讲授内容：龙格库塔方法、收敛性与稳定性 三、教学重点难点 1．教学重点：龙格库塔方法、收敛性与稳定性。 2. 教学难点：收敛性与稳定性。 四、教学中应注意的问题 多媒体课堂教学为主。适当提问，加深学生对概念的理解。 五、正文 龙格－库塔方法 引言 龙格－库塔方法的基本思想 '11()() ()()()(,())n n n n y x y x y y x y x hf y h ξξξ++-=?=+ 令*(,())K f y ξξ=称为区间1[,]n n x x +上的平均斜率，只要对*K 提供一种算法，即可推导出一种计算公式。 欧拉公式只是取点n x 的斜率作为区间1[,]n n x x +的平均斜率*K ，精度自然很低。 考察改进的欧拉公式 11211 ()22 n n y y h k k +=++1(,)n n k f x y =21(,)n n k f x h y hk =++ 它利用n x ，1n x +两点的斜率取算术平均，1n x +处斜率通过已知信息n y 用欧拉公式预报得到。

可以考虑设法在1[,]n n x x +多取几个点的斜率值，将它们加权平均作为区间 1[,]n n x x +的平均斜率*K 。这就是龙格－库塔方法的基本思想。 1、二阶龙格－库塔方法 考察1[,]n n x x +内一点,01n p n x x hp p +=+<≤，希望通过,n n p x x +两个点的斜率 12,K K 加权平均得到*K ，即令112((1))n n y y h K K λλ+=+-+ 取1(,)n n K f x y =，如何预报n p x +处斜率2K ？ 仿照改进的欧拉公式，先用欧拉公式预测()n p y x +的值n p y +： 1n p n y y phK +=+ 然后用n p y +计算2(,)n p n p K f x y ++=，从而得 112121[(1)](,)(,) n n n n n p n y y h K K K f x y K f x y phK λλ++=+-+==+ 适当选取,p λ，使上述公式具有较高得精度。假定()n n y y x =，分别将12 ,K K 泰勒展开： '1(,)()n n n K f x y y x == 212 ' '' 2 (,) (,)[(,)(,)(,)]()()()() n p n n n x n n n n y n n n n K f x y pkK f x y ph f x y f x y f x y O h y x phy x O h +=+=+++=++ 代入得 '2''31()()()()n n n n y y x hy x ph y x O h λ+=+++ 按泰勒展开2' '' 31()()()()2 n n n n h y y x hy x y x O h +=+++ 比较得，只要1 2 p λ= ，公式截断误差为3()O h 特别，当1,1/2p λ==，就是改进的欧拉公式， 改取1,1/2p λ==， 12n n y y hk +=+，1(,)n n k f x y =，21(,)22 n n h h k f x y k =++

四阶龙格库塔法的编程(赵)

例题一 四阶龙格-库塔法的具体形式为： 1.1程序: /*e.g: y'=t-y,t∈[0,1] /*y(0)=0 /*使用经典四阶龙格-库塔算法进行高精度求解 /* y(i＋1)＝yi＋( K1+ 2*K2 +2*K3+ K4)/6 /* K1=h*f(ti,yi) /* K2=h*f(ti+h/2,yi+K1*h/2) /* K3=h*f(ti+h/2,yi+K2*h/2) /* K4=h*f(ti+h,yi+K3*h) */ #include #include #define N 20 float f(float d,float p) //要求解的微分方程的右部的函数e.g: t-y { float derivative; derivative=d-p; return(derivative); } void main() { float t0; //范围上限

float t; //范围下限 float h; //步长 float nn; //计算出的点的个数,即迭代步数 int n; //对nn取整 float k1,k2,k3,k4; float y[N]; //用于存放计算出的常微分方程数值解 float yy; //精确值 float error;//误差 int i=0,j; //以下为函数的接口 printf("input t0:"); scanf("%f",&t0); printf("input t:"); scanf("%f",&t); printf("input y[0]:"); scanf("%f",&y[0]); printf("input h:"); scanf("%f",&h); // 以下为核心程序 nn=(t-t0)/h; printf("nn=%f\n",nn); n=(int)nn; printf("n=%d\n",n); for(j=0;j<=n;j++) { yy=t0-1+exp(-t0); //解析解表达式 error=y[i]-yy; //误差计算 printf("y[%f]=%f yy=%f error=%f\n",t0,y[i],yy,error);//结果输出k1=h*f(t0,y[i]); //求K1 k2=h*f((t0+h/2),(y[i]+k1*h/2)); //求K2 k3=h*f((t0+h/2),(y[i]+k2*h/2)); //求K3

四阶龙格库塔法原理C代码

/** ***四阶Runge-Kutta法*** 经典格式： y(n+1) = y(n) + h/6 ( K1 + 2*K2 + 2*K3 + K4 ) K1 = f( x(n) , y(n) ) K2 = f( x(n+1/2) , y(n) + h/2*K1 ) K3 = f( x(n+1/2) , y(n) + h/2*K2 ) K4 = f( x(n+1) , y(n) + h*K3 ) Runge-Kutta法是基于泰勒展开方法，因而需要所求解具有较好的光滑性。 属性：差分方法 《数值分析简明教程》-2 Editon -高等教育出版社- page 105 算法流程图代码维护：2005.6.14 DragonLord **/ #include #include #include /* 举例方程： y'= y - 2*x / y ( 0>x0>>y0>>h>>N) { int n=0;

for(;n

龙格库塔方法matlab实现

龙格库塔方法matlab实现~ function ff=rk(yy,x0,y0,h,a,b)%yy为y的导函数,x0,y0,为初值,h为步长,a,b为区间 c=(b-a)/h+1;i1=1; %c为迭代步数；i1为迭代步数累加值 y=y0;z=zeros(c,6); %z生成c行，5列的零矩阵存放结果； %每行存放c次迭代结果，每列分别存放k1~k4及y的结果 for x=a:h:b if i1<=c k1=feval(yy,x,y); k2=feval(yy,x+h/2,y+(h*k1)/2); k3=feval(yy,x+h/2,y+(h*k2)/2); k4=feval(yy,x+h,y+h*k3); y=y+(h/6)*(k1+2*k2+2*k3+k4); z(i1,1)=x;z(i1,2)=k1;z(i1,3)=k2;z(i1,4)=k3;z(i1,5)=k4;z(i1,6)=y; i1=i1+1; end end fprintf(‘结果矩阵，第一列为x(n),第二列~第五列为k1~k4，第六列为y（n+1）的结果') z %在命令框输入下列语句 %yy=inline('x+y'); %>> rk(yy,0,1,0.2,0,1) %将得到结果 %结果矩阵，第一列为x(n),第二列~第五列为k1~k4第六列为y（n+1）的结果 %z = % 0 1.0000 1.2000 1.2200 1.4440 1.2428 % 0.2000 1.4428 1.6871 1.7115 1.9851 1.5836 % 0.4000 1.9836 2.2820 2.3118 2.6460 2.0442 % 0.6000 2.6442 3.0086 3.0451 3.4532 2.6510 % 0.8000 3.4510 3.8961 3.9407 4.4392 3.4365 % 1.0000 4.4365 4.9802 5.0345 5.6434 4.4401

标准四阶龙格——库塔法

实验名：常微分方程数值解法 实习目的： （1） 通过实习进一步掌握标准四阶龙格——库塔法的基本思想； （2） 通过对标准四阶龙格——库塔法的调试练习，进一步体会其特点； （3） 通过实习进一步掌握标准四阶龙格——库塔法的计算步骤，并能灵活应用； （4） 通过上机调试运行，逐步培养解决实际问题的编程能力。 实习要求： （1） 熟悉Turbo C 的编译环境； （2） 实习前复习标准四阶龙格——库塔法的基本思想和过程； （3） 实习前复习标准四阶龙格——库塔法的计算步骤。 实习设备： （1） 硬件设备：单机或网络环境下的微型计算机一台; （2） 软件设备：DOS3.3以上操作系统，Turbo C2.0编译器。 实习内容： 标准四阶龙格——库塔法： （1）使用标准四阶龙格——库塔法求解初值问题 的数值求解。 （2）要求： 请写出程序的运行结果： 程序代码： #include "stdio.h" #include "conio.h" float func(float x,float y) { return(2*x*y); } float runge_kutta(float x0,float xn,float y0,int n) { float x,y,y1,y2,h,xh; float d1,d2,d3,d4; int i; x=x0; y=y0; h=(xn-x0)/n; for(i=1;i<=n;i++) 1(0)y 1x 0 2xy y =≤≤='

{ xh=x+h/2; d1=func(x,y); d2=func(xh,y+h*d1/2.0); d3=func(xh,y+h*d2/2.0); d4=func(xh,y+h*d3); y=y+h*(d1+2*d2+2*d3+d4)/6.0; x=x0+i*h; } return(y); } void main() { float x0,xn,y0,e; int n; printf("\ninput n:\n"); scanf("%d",&n); printf("input x0,xn:\n"); scanf("%f%f",&x0,&xn); printf("input y0:\n"); scanf("%f",&y0); e=runge_kutta(x0,xn,y0,n); printf("y(%f)=%6.6f",y0,e); } 运行结果: （3）思考题： 标准四阶龙格——库塔法的基本思想是什么? 龙格和库塔提出了一种间接地运用Taylor公式的方法，即利用y(x)在若干个待定点上的函数值和导数值做出线性组合式，选取适当系数使这个组合式进Taylor展开后与y(xi+1)的Taylor 展开式有较多的项达到一致，从而得出较高阶的数值公式，这就是龙格—库塔法的基本思想。

龙格库塔方法及其matlab实现

龙格-库塔方法及其matlab实现 摘要：本文的目的数值求解微分方程精确解，通过龙格-库塔法，加以利用matlab为工具 达到求解目的。龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法，用于数值求解微分方程。MatLab软件是由美国Mathworks公司推出的用于数值计算和图形 处理的科学计算系统环境。MatLab是英文MATrix LABoratory（矩阵实验室）的缩写。在MratLab环境下，用户可以集成地进行程序设计、数值计算、图形绘制、输入输出、文件 管理等各项操作。 关键词：龙格-库塔 matlab 微分方程 1.前言 1.1：知识背景 龙格－库塔法（Runge-Kutta）是用于非线性常微分方程的解的重要的一类隐式或显式迭代法。这些技术由数学家卡尔·龙格和马丁·威尔海姆·库塔于1900年左右发明。通常所说的龙格库塔方法是相对四阶龙格库塔而言的，成为经典四阶龙格库塔法。该方法具有精度高，收敛，稳定，计算过程中可以改变步长不需要计算高阶导数等优点，但是仍需计算在 一些点上的值，比如四阶龙格-库塔法没计算一步需要计算四步，在实际运用中是有一定复杂性的。 Matlab是在20世纪七十年代后期的事：时任美国新墨西哥大学计算机科学系主任的Cleve Moler教授出于减轻学生编程负担的动机，为学生设计了一组调用LINPACK和EISPACK库 程序的“通俗易用”的接口，此即用FORTRAN编写的萌芽状态的MATLAB。 经几年的校际流传，在Little的推动下，由Little、Moler、Steve Bangert合作，于1984年成立了MathWorks公司，并把MATLAB正式推向市场。从这时起，MATLAB的内核 采用C语言编写，而且除原有的数值计算能力外，还新增了数据图视功能。 MATLAB以商品形式出现后，仅短短几年，就以其良好的开放性和运行的可靠性， 使原先控制领域里的封闭式软件包（如英国的UMIST，瑞典的LUND和SIMNON，德国的KEDDC）纷纷淘汰，而改以MATLAB为平台加以重建。在时间进入20世纪九十年代的时候，MATLAB已经成为国际控制界公认的标准计算软件。 到九十年代初期，在国际上30几个数学类科技应用软件中，MATLAB在数值计算方面独占 鳌头，而Mathematica和Maple则分居符号计算软件的前两名。Mathcad因其提供计算、 图形、文字处理的统一环境而深受中学生欢迎。 1.2研究的意义 精确求解数值微分方程，对龙格库塔的深入了解与正确运用，主要是在已知方程导数和初 值信息，利用计算机仿真时应用，省去求解微分方程的复杂过程。利用matlab强大的数值计算功能，省去认为计算的过程，达到快速精确求解数值微分方程。在实际生活中可以利 用龙格库塔方法和matlab的完美配合解决问题。 1.3研究的方法 对实例的研究对比，实现精度的要求，龙格库塔是并不是一个固定的公式，所以只是对典 型进行分析

龙格库塔法例题

四阶龙格一库塔法 通常所说的龙格一库塔法是指四阶而言的．我们可以仿二阶、三阶的情形推导出常用的标准四阶龙格一库塔公式 (9.22) 公式(9.22)的局部截断误差的阶为． 龙格一库塔法具有精度高，收敛，稳定(在一定的条件下)，计算过程中可以改变步长，不需要计算高阶导数值等优点．但仍需计算在一些点上的值，如四阶龙格－库塔法每计算一步需要算四次 的值，这给实际计算带来一定的复杂性，因此，多用来计算“表头”．(即开始几点的近似值)．例3.用标准龙格一库塔法求初值问题 在处的解． 解因与．若应用标准龙格一库塔方法公式(9.22)计算，对于n=0时，则有

于是得 这个值与准确解在处的值已十分接近．再对n=1，2，3，4应用式(9.22)计算，具体计算结果如表3所示：

例3写出用四阶龙格――库塔法求解初值问题 的计算公式，取步长h＝0.2计算y(0.4)的近似值。至少保留四位小数。 解此处f(x，y)＝8－3y，四阶龙格――库塔法公式为 其中κ1＝f(x k，y k)；κ2＝f(x k+0.5h，y k+0.5hκ1)； κ3＝f(x k+0.5h，y k+0.5hκ2)；κ4＝f(x k+h，y k+hκ3) 本例计算公式为： 其中κ1＝8－3y k；κ2＝5.6－2.1y k； κ3＝6.32－2.37y k；κ4＝4.208－1.578y k ＝1.2016＋0.5494y k (k＝0，1，2，…) 当x0＝0，y0＝2， y(0.2)≈y1＝1.2016＋0.5494y0＝1.2016＋0.5494×2＝2.3004 y(0.4)≈y2＝1.2016＋0.5494y1＝1.2016＋0.5494×2.3004＝2.4654

常微分方程组的四阶RungeKutta龙格库塔法matlab实现

常微分方程组的四阶Runge-Kutta方法1.问题： 1.1若用普通方法-----仅适用于两个方程组成的方程组 编程实现： 创建M 文件： function R = rk4(f,g,a,b,xa,ya,N) %UNTITLED2 Summary of this function goes here % Detailed explanation goes here %x'=f(t,x,y) y'=g(t,x,y) %N为迭代次数 %h为步长 %ya,xa为初值 f=@(t,x,y)(2*x-0.02*x*y);

g=@(t,x,y)(0.0002*x*y-0.8*y); h=(b-a)/N; T=zeros(1,N+1); X=zeros(1,N+1); Y=zeros(1,N+1); T=a:h:b; X(1)=xa; Y(1)=ya; for j=1:N f1=feval(f,T(j),X(j),Y(j)); g1=feval(g,T(j),X(j),Y(j)); f2=feval(f,T(j)+h/2,X(j)+h/2*f1,Y(j)+g1/2); g2=feval(g,T(j)+h/2,X(j)+h/2*f1,Y(j)+h/2*g1); f3=feval(f,T(j)+h/2,X(j)+h/2*f2,Y(j)+h*g2/2); g3=feval(g,T(j)+h/2,X(j)+h/2*f2,Y(j)+h/2*g2); f4=feval(f,T(j)+h,X(j)+h*f3,Y(j)+h*g3); g4=feval(g,T(j)+h,X(j)+h*f3,Y(j)+h*g3); X(j+1)=X(j)+h*(f1+2*f2+2*f3+f4)/6; Y(j+1)=Y(j)+h*(g1+2*g2+2*g3+g4)/6; R=[T' X' Y']; end 情况一：对于x0=3000，y0=120 控制台中输入： >> rk4('f','g',0,10,3000,120,10) 运行结果： ans = 1.0e+003 * 0 3.0000 0.1200 0.0010 2.6637 0.0926 0.0020 3.7120 0.0774 0.0030 5.5033 0.0886 0.0040 4.9866 0.1193 0.0050 3.1930 0.1195 0.0060 2.7665 0.0951 0.0070 3.6543 0.0799 0.0080 5.2582 0.0884 0.0090 4.9942 0.1157 0.0100 3.3541 0.1185 数据:

龙格-库塔法MATLAB

1. matlab 新建.m文件，编写龙格-库塔法求解函数 function [x,y]=runge_kutta1(ufunc,y0,h,a,b)%参数表顺序依次是微分方程组的函数名称，初始值向量，步长，时间起点，时间终点（参数形式参考了ode45函数） n=floor((b-a)/h); %求步数 x(1)=a;%时间起点 y(:,1)=y0;%赋初值，可以是向量，但是要注意维数 for ii=1:n x(ii+1)=x(ii)+h; k1=ufunc(x(ii),y(:,ii)); k2=ufunc(x(ii)+h/2,y(:,ii)+h*k1/2); k3=ufunc(x(ii)+h/2,y(:,ii)+h*k2/2); k4=ufunc(x(ii)+h,y(:,ii)+h*k3); y(:,ii+1)=y(:,ii)+h*(k1+2*k2+2*k3+k4)/6; %按照龙格库塔方法进行数值求解 end 2.另外再新建一个.,m文件，定义要求解的常微分方程函数 function dx=fun1(t,x) dx =zeros(2,1);%初始化列向量 dx(1) =0.08574*x(2)-1.8874*x(1)-20.17; dx(2) =1.8874*x(1)-0.08574*x(2)+115.16; 3,再新建一个.m文件，利用龙格-库塔方法求解常微分方程 [T1,F1]=runge_kutta1(@fun1,[46.30 1296 ],1,0,20); %求解步骤2定义的fun1常微分方程，@fun1是调用其函数指针，从0到20，间隔为1 subplot(122) plot(T1,F1)%自编的龙格库塔函数效果 title('自编的龙格库塔函数') grid 运行步骤3文件即可得到结果，F1为估测值 或者可以调用matlab自带函数ode45求解 方法如下：

1、经典四阶龙格库塔法解一阶微分方程组

陕西科技大学 数值计算课程设计任务书 理学院信息与计算科学/应用数学专业信息08/数学08 班级学生： 题目：典型数值算法的C++语言程序设计 课程设计从2010 年 5 月17日起到2010 年 6 月18 日 1、课程设计的内容和要求（包括原始数据、技术要求、工作要求等）： 每人需作10个算法的程序、必做6题、自选4题。 对每个算法要求用C++语言进行编程。 必选题： 1、经典四阶龙格库塔法解一阶微分方程组 2、高斯列主元法解线性方程组 3、牛顿法解非线性方程组 4、龙贝格求积分算法 5、三次样条插值算法（压紧样条）用C++语言进行编程计算 依据计算结果，用Matlab画图并观察三次样条插值效果。 6、M次多项式曲线拟合，据计算结果，用Matlab画图并观察拟合效果。 自选题：自选4道其他数值算法题目.每道题目重选次数不得超过5次. 2、对课程设计成果的要求〔包括图表、实物等硬件要求〕： 1）提交课程设计报告 按照算法要求，用C++语言设计和开发应用程序，提交由算法说明；程序设计说明；系统技术文档(包括系统各模块主要流程图，软件测试方案与测试记录、软件调试和修改记录、测试结论、运行情况记录)，系统使用说明书，源程序代码为附录构成的课程设计报告。 2）课程设计报告版式要求

打印版面要求：A4纸，页边距：上2cm，下2cm，左2.5cm、右2cm；字体：正文宋体、小四号；行距：固定值20；页眉1.5cm ，页脚1.75cm；页码位于页脚居中打印；奇数页页眉“数值计算课程设计”，偶数页页眉“算法名称”，页眉宋体小5号；段落及层次要求：每节标题以四号黑体左起打印（段前段后各0.5行），节下为小节，以小四号黑体左起打印（段前段后各0.5行）。换行后以小四号宋体打印正文。节、小节分别以1、1.1、1.1.1依次标出，空一字符后接各部分的标题。 当论文结构复杂，小节以下的标题，左起顶格书写，编号依次用（1）、（2）……或1）、2）……顺序表示。字体为小四号宋体。 对条文内容采用分行并叙时，其编号用（a）、（b）……或a）、b）……顺序表示，如果编号及其后内容新起一个段落，则编号前空两个中文字符。3）设计报告装订顺序与规范 封面 数值计算课程设计任务书 目录 数值计算设计课程设计报告正文 设计体会及今后的改进意见 参考文献（资料） 左边缘装订 3 指导教师：日期： 教研室主任：日期：

二阶龙格库塔方法

2012-2013（1）专业课程实践论文二阶Runge-Kutta方法 董文峰，0818180123，R数学08-1班

由改进的Euler 方法得到： ()) ,(),(21121 211 K y h x hf K y x hf K K K y y i i i i i i ++==++=?????+ 凡满足条件式有一簇形如上式的计算格式，这些格式统称为二阶龙格—库塔格式。因此改进的欧拉格式是众多的二阶龙格—库塔法中的一种特殊格式。 若取1,0,2121212=== =c c b a ，就是另一种形式的二阶龙格 - 库塔公式。 ??????? ++==+=+)2,2() ,(12121K h y h x f K y x f K hK y y n n n n n n （1） 此计算公式称为变形的二阶龙格—库塔法。 二级龙格-库塔方法是显式单步式,每前进一步需要计算两个函数值。由上面的讨论可知,适当选择四个参数y0,a,b,n ，可使每步计算两次函数值的二阶龙格-库塔方法达到二阶精度。下面以式子（1）为依据利用VC++6.0编译程序进行问题的求解。

#include #include /*n表示几等分，n+1表示他输出的个数*/ int RungeKutta(double y0,double a,double b,int n,double *x,double *y,double (*function)(double,double)) { double h=(b-a)/n,k1,k2; int i; x[0]=a; y[0]=y0; for(i=0;i

实验五 变步长的龙哥库塔法

实验五变步长的龙格库塔法 一、实验目的 1、了解变步长的龙哥库塔法的特点及具体实现过程。 2、运用vc6.0软件编写龙哥库塔法的算法程序。 3、采用编制的程序，对实际问题进行数值模拟计算。 4、将数值计算结果与精确解进行对比分析，讨论计算误差。 二、实验内容 1、根据变步长的龙哥库塔法法算法的特点，设计程序的流程 本实验采用经典的变步长的龙哥库塔法法格式（四阶变步长的龙哥库塔法法格式） 2、用vc6.0语言编写变步长的龙哥库塔法法算法程序 在编写程序时，充分考虑到程序的交互性和实用性。本程序可以实现计算机主动提示操作步骤和输出结果； 输入待计算的信息后，计算机自动生成图像，并且将得到的结果也显示在图像中，便于直观观察； 3、算法的源程序 #include "stdio.h" #include "math.h" double f(double,double); main() { double y0=1,k1,k2,k3,k4,x0=0,h=0.1,h1=0.05,q,tol=0.1,y00,y01; int i,n=10; k1=f(x0,y0); k2=f(x0+h/2,y0+h/2*k1); k3=f(x0+h/2,y0+h/2*k2); k4=f(x0+h/2,y0+h*k3); y0=y0+h/6*(k1+2*k2+2*k3+k4); y00=y0; printf("the y00 is: %lf\n",y00); /*步长为h时从x0出发求一步得y00*/ k1=f(x0,y0); k2=f(x0+h1/2,y0+h1/2*k1); k3=f(x0+h1/2,y0+h1/2*k2); k4=f(x0+h1/2,y0+h1*k3); y0=y0+h1/6*(k1+2*k2+2*k3+k4); x0=x0+h1; k1=f(x0,y0); k2=f(x0+h1/2,y0+h1/2*k1); k3=f(x0+h1/2,y0+h1/2*k2);

四阶龙格——库塔法

2013-2014（1）专业课程实践论文题目：四阶龙格—库塔法

一、算法理论 由定义可知，一种数值方法的精度与局部截断误差()p o h 有关，用一阶泰勒展开式近似函数得到欧拉方法，其局部截断误差为一阶泰勒余项2()o h ，故是一阶方法，完全类似地若用p 阶泰勒展开式 2 '''() 11()()()......()() 2!!p p p n n n n n h h y y x hy x y x y x O h p ++=+++++ 进行离散化，所得计算公式必为p 阶方法，式中 '''''()(,),()(,)(,)(,).... x y x f x y y x f x y f x y f x y ==++ 由此，我们能够想到，通过提高泰勒展开式的阶数，可以得到高精度的数值方法，从理论上讲，只要微分方程的解()y x 充分光滑，泰勒展开方法可以构造任意的有限阶的计算公式，但事实上，具体构造这种公式往往相当困难，因为符合函数(,())f x y x 的高阶导数常常是很烦琐的，因此，泰勒展开方法一般不直接使用，但是我们可以间接使用泰勒展开方法，求得高精度的计算方法。 首先，我们对欧拉公式和改进欧拉公式的形式作进一步的分析。 如果将欧拉公式和改进的欧拉公式改写成如下的形式： 欧拉公式 {111(,)n n n n y y hK K f x y +==+ 改进的欧拉公式 11211()22 n n y y h K K +=++， 1(,)n n K f x y =， 21(,)n n K f x h y hK =++。 这两组公式都是用函数(,)f x y 在某些点上的值的线性组合来计算1()n y x +的近似值1n y +，欧拉公式每前进一步，就计算一次(,)f x y 的值。另一方面它是1()n y x +在n x 处的一阶泰勒展开式，因而是一阶方法。改进的欧拉公式每前进一步，

龙格库塔法-原理及程序实现

龙格库塔法一、基本原理：

可以得出四阶龙格－库塔公式，也就是在工程中应用广泛的经典龙格－库塔算法，即： yi＋１＝yi＋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) 通常所说的龙格-库塔法就是指四阶——龙格库塔法，我们可以仿二阶、三阶的情形推导出常用的标准四阶龙格-库塔法公式。 （1） 计算公式（1）的局部截断误差是。 龙格-库塔法具有精度高，收敛，稳定（在一定条件下），计算过程中可以改变步长，不需要计算高阶导数等优点，但仍需计算 在一些点上的值，如四阶龙格-库塔法每计算一步需要计算四次 的值，这给实际计算带来一定的复杂性，因此，多用来计算“表头”。 二、小程序 #include #include

#define f(x,y) (-1*(x)*(y)*(y)) void main(void) { double a,b,x0,y0,k1,k2,k3,k4,h; int n,i; printf("input a,b,x0,y0,n:"); scanf("%lf%lf%lf%lf%d",&a,&b,&x0,&y0,&n); printf("x0\ty0\tk1\tk2\tk3\tk4\n"); for(h=(b-a)/n,i=0;i!=n;i++) { k1=f(x0,y0); k2=f(x0+h/2,y0+k1*h/2); k3=f(x0+h/2,y0+k2*h/2); k4=f(x0+h,y0+h*k3); printf("%lf\t%lf\t",x0,y0); printf("%lf\t%lf\t",k1,k2); printf("%lf\t%lf\n",k3,k4); y0+=h*(k1+2*k2+2*k3+k4)/6; x0+=h; } printf("xn=%lf\tyn=%lf\n",x0,y0); }

四阶龙格库塔法解微分方程

四阶龙格库塔法解微分方程 一、四阶龙格库塔法解一阶微分方程 ①一阶微分方程： cos y t & ，初始值y(0)=0，求解区间[0 10]。 MATLAB 程序： %%%%%%%%%%% 四阶龙哥库塔法解一阶微分方程 %%%%%%%%%%% y'=cost %%%%%%%%%%% y(0)=0, 0≤t ≤10，h=0.01 %%%%%%%%%%% y=sint h=0.01; hf=10; t=0:h:hf; y=zeros(1,length(t)); y(1)=0; F=@(t,y)(cos(t)); for i=1:(length(t)-1) k1=F(t(i),y(i)); k2=F(t(i)+h/2,y(i)+k1*h/2); k3=F(t(i)+h/2,y(i)+k2*h/2); k4=F(t(i)+h,y(i)+k3*h); y(i+1)=y(i)+1/6*(k1+2*k2+2*k3+k4)*h;

end subplot(211) plot(t,y,'-') xlabel('t'); ylabel('y'); title('Approximation'); span=[0,10]; p=y(1); [t1,y1]=ode45(F,span,p); subplot(212) plot(t1,y1) xlabel('t'); ylabel('y'); title('Exact');

图1 ②一阶微分方程： ()22*/x t x x t =-&& ，初始值x(1)=2，求解区间[1 3]。 MATLAB 程序： %%%%%%%%%%% 四阶龙哥库塔法解微分方程 %%%%%%%%%%% x'(t)=(t*x-x^2)/t^2 %%%%%%%%%%% x(1)=2, 1≤t ≤3, h=1/128 %%%%%%%%%%% 精确解：x(t)=t/(0.5+lnt) h=1/128; %%%%% 步长 tf=3; t=1:h:tf;

四阶龙格库塔法解一阶二元微分方程

四阶龙格库塔法解一阶二元微分方程 //dxi/dt=c*(xi-xi^3/3+yi)+K*(X-xi)+c*zi //dyi/dt=(xi-b*yi+a)/c //i=1,2,3 //X=sum(xi)/N #include #include #include #include #define N 1000 //定义运算步数； #define h 0.01 //定义步长； float a,b,c;//定义全局变量常数a，b，c //定义微分方程： double fx(double x[],double dx,double y[],double dy,double z[],int i,double k,double xavg){ int j; double xi,yi; xi=x[i]+dx; yi=y[i]+dy; return c*(xi-pow(xi,3)/3+yi)+k*(xavg-xi)+c*z[i]; } double fy(double x[],double dx,double y[],double dy,int i){ double xi,yi; xi=x[i]+dx; yi=y[i]+dy; return (xi-b*yi+a)/c; } void main(){ double Kx[3][4],Ky[3][4],x[3]={1,2,3},y[3]={2,3,4},xavg,k=0;//定义x,y的初值；double z[3]={0}; int i,j,m,n,S; FILE *fp1,*fp; fp=fopen("sjy.txt","w"); fp1=fopen("sjykxy.txt","w"); fprintf(fp1,"k\tx1\tx2\tx3\ty1\ty2\ty3\n"); if(fp==NULL||fp1==NULL){ printf("Failed to open file.\n"); getch(); return;