数值计算方法第4次实验--龙格库塔公式

1.完成复合梯形、复合辛普森求积公式，用变步长、事后误差估计的方法完成P145例题6.3.1对应表6-3

的计算，课后作业题7.

%复合梯形求积公式的MATLAB实现%姓名：王定

%学号：1306034248

%中北大学仪器与电子学院

function T_n=ComTrapezoidal(a,b,n)

%a为积分上限

%b为积分下限

%n为划分区间的个数

h=(b-a)/n;

for(k=0:n)

x(k+1)=a+k*h;

if(x(k+1)==0)

x(k+1)=10^(-10);

end

end

I_1=h/2*(f(x(1))+f(x(n+1)));

for(i=2:n)

F(i)=h*f(x(i));

end

I_2=sum(F);

I_n=I_1+I_2

%****************************% %复合辛普森求积公式的MATLAB实现

function I=ComSimpson(a,b,n)

n=2;

h=(b-a)/2;

I1=0;

I2=(f(a)+f(b))/h;

eps=1.0e-4;

while(abs(I2-I1)>eps)

n=n+1;

h=(b-a)/n;

I1=I2;

I2=0;

for(i=0:n-1)

x=a+h*i;

x1=x+h;

I2=I2+h/6*(f(x)+4*f((x+x1)/2)+f(x1));

end

end

I=I2 %变步长梯形求积法的MATLAB实

现

function AdaptiveTrapezoidal(a,b,eps)

m=1

t0=0;

t2=(f(a)+f(b))/4+f((a+b)/2)

while(abs(t2-t0)>eps)

m=m+1

p=t2;

t1=0;

n=2^m;

h=(b-a)/n;

for(k=0:(n-1))

t1=t1+h*((f(a+k*h)+f(a+(k+1)*h))/4+

f(a+(k+1/2)*h)/2);

end

t2=t1

t0=p;

end

I=t2

%**************************%

>>

AdaptiveTrapezoidal(1,2,0.000001)

m = 1

t2 = 3.03948481584447

m = 2

t2 = 2.02304986763725

m = 3

t2 = 2.02080858246806

m = 4

t2 = 2.02024624995310

m = 5

t2 = 2.02010553934348

m = 6

t2 = 2.02007035369492

m = 7

t2 = 2.02006155678257

m = 8

t2 = 2.02005935752321

m = 9

t2 = 2.02005880770641

I = 2.02005880770641

%课后作业题7.

>>

AdaptiveTrapezoidal(1,

3,0.0001)

m = 1

t2 =

7.99930630234471

m = 2

t2 =

10.84204346755743

m = 3

t2 =

10.92309388961378

m = 4

t2 =

10.94339842118680

m = 5

t2 =

10.94847716722413

m = 6

t2 =

10.94974701694155

m = 7

t2 =

10.95006448956964

m = 8

t2 =

10.95014385836406

I =

10.95014385836406

龙格库塔方法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

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

龙格库塔积分算法

龙格库塔法 龙格库塔法是常用于模拟常微分方程的解的重要的一类隐式或显式迭代法。这些技术由数学家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）

Matlab中龙格-库塔(Runge-Kutta)方法原理及实现

函数功能编辑本段回目录 ode是专门用于解微分方程的功能函数，他有ode23,ode45,ode23s等等，采用的是Runge-Kutta算法。ode45表示采用四阶，五阶runge-kutta单步算法,截断误差为(Δx)3。解决的是Nonstiff(非刚性)的常微分方程.是解决数值解问题的首选方法,若长时间没结果,应该就是刚性的,换用ode23来解. 使用方法编辑本段回目录 [T,Y] = ode45(odefun,tspan,y0) odefun 是函数句柄,可以是函数文件名,匿名函数句柄或内联函数名 tspan 是区间[t0 tf] 或者一系列散点[t0,t1,...,tf] y0 是初始值向量 T 返回列向量的时间点 Y 返回对应T的求解列向量 [T,Y] = ode45(odefun,tspan,y0,options) options 是求解参数设置,可以用odeset在计算前设定误差,输出参数,事件等 [T,Y,TE,YE,IE] =ode45(odefun,tspan,y0,options) 在设置了事件参数后的对应输出 TE 事件发生时间 YE 事件解决时间 IE 事件消失时间 sol =ode45(odefun,[t0 tf],y0...) sol 结构体输出结果 应用举例编辑本段回目录 1 求解一阶常微分方程

程序: 一阶常微分方程 odefun=@(t,y) (y+3*t)/t^2; %定义函数 tspan=[1 4]; %求解区间 y0=-2; %初值 [t,y]=ode45(odefun,tspan,y0); plot(t,y) %作图 title('t^2y''=y+3t,y(1)=-2,1

Matlab中龙格-库塔(Runge-Kutta)方法原理及实现

Matlab中龙格-库塔(Runge-Kutta)方法原理及实现 龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。该算法是构建在数学支持的基础之上的。龙格库塔方法的理论基础来源于泰勒公式和使用斜率近似表达微分，它在积分区间多预计算出几个点的斜率，然后进行加权平均，用做下一点的依据，从而构造出了精度更高的数值积分计算方法。如果预先求两个点的斜率就是二阶龙格库塔法，如果预先取四个点就是四阶龙格库塔法。一阶常微分方程可以写作：y'=f(x,y),使用差分概念。 (Yn+1-Yn)/h= f(Xn,Yn)推出（近似等于，极限为Yn'） Yn+1=Yn+h*f(Xn,Yn) 另外根据微分中值定理，存在0

所以，为了更好更准确地把握时间关系，应自己在理解龙格库塔原理的基础上，编写定步长的龙格库塔函数，经过学习其原理，已经完成了一维的龙格库塔函数。 仔细思考之后，发现其实如果是需要解多个微分方程组，可以想象成多个微分方程并行进行求解，时间，步长都是共同的，首先把预定的初始值给每个微分方程的第一步，然后每走一步，对多个微分方程共同求解。想通之后发现，整个过程其实很直观，只是不停的逼近计算罢了。编写的定步长的龙格库塔计算函数： 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; %按照龙格库塔方法进行数值求解

龙格库塔法求微方程matlab

龙格—库塔方法求解微分方程初值问题 （数学1201+41262022+陈晓云） 初值问题： y x x -+=2dx dy ，10≤≤x 1)0(y = 四阶龙格-库塔公式： ()y x K n n ,f 1= ????? ? ??+=+K h y x K n h n 122f ,2 ??? ??++=K y x f K h n h n 232,2 ()K h y h x f K n n 34,++= ()K K K K y y h n 4 3211n 226++++=+ 程序： 1）建立四阶龙格-库塔函数 function [ x,y ] = nark4( dyfun,xspan,y0,h ) % dyfun 为一阶微分方程的函数；y0为初始条件；xspan 表示x 的区间；h 为区间的步长； x=xspan(1):h:xspan(2); y(1)=y0; for n=1:length(x)-1 k1=feval(dyfun,x(n),y(n)); k2=feval(dyfun,x(n)+h/2,y(n)+h/2*k1); k3=feval(dyfun,x(n)+h/2,y(n)+h/2*k2); k4=feval(dyfun,x(n+1),y(n)+h*k3); y(n+1)=y(n)+h*(k1+k2*2+2*k3+k4)/6; end x=x;y=y;

2）执行程序（m文件） dyfun=inline('x^2+x-y'); [x,y1]=nark4(dyfun,[0,1],1,0.1); x=0:0.1:1; Format long y2=x.^2-x+1 R4=y2-y1 [x',y1',y2',R4'] y2=dsolve('Dy=x^2+x-y','y(0)=1','x') plot(x,y1,'b*-') hold on y3=inline('x^2-x+1') fplot(y3,[0,1],'ro-') legend('R-K4','解析解') 3）执行结果 ans = X RK4近似值解析值 0 1.000000000000000 1.000000000000000 0.100000000000000 0.910000208333333 0.910000000000000 0.200000000000000 0.840000396841146 0.840000000000000 0.300000000000000 0.790000567410084 0.790000000000000 0.400000000000000 0.760000721747255 0.760000000000000 0.500000000000000 0.750000861397315 0.750000000000000 0.600000000000000 0.760000987757926 0.760000000000000 0.700000000000000 0.790001102093746 0.790000000000000 0.800000000000000 0.840001205549083 0.840000000000000 0.900000000000000 0.910001299159352 0.910000000000000 1.000000000000000 1.000001383861433 1.000000000000000

龙格库塔法例题

四阶龙格一库塔法 通常所说的龙格一库塔法是指四阶而言的．我们可以仿二阶、三阶的情形推导出常用的标准四阶龙格一库塔公式 (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

龙格库塔方法及其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研究的方法 对实例的研究对比，实现精度的要求，龙格库塔是并不是一个固定的公式，所以只是对典 型进行分析

二阶龙格库塔方法

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

欧拉法与龙格库塔法比较分析

解微分方程的欧拉法，龙格-库塔法简单实例比较 欧拉方法（Euler method)用以对给定初值的常微分方程(即初值问题)求解分为前EULER 法、后退EULER 法、改进的EULER 法。 缺点： 欧拉法简单地取切线的端点作为下一步的起点进行计算，当步数增多时，误差会因积累而越来越大。因此欧拉格式一般不用于实际计算。 改进欧拉格式（向前欧拉公式）： 为提高精度，需要在欧拉格式的基础上进行改进。采用区间两端的斜率的平均值作为直线方程的斜率。改进欧拉法的精度为二阶。 算法： 微分方程的本质特征是方程中含有导数项，数值解法的第一步就是设法消除其导数值。对于常微分方程： (,)dy f x y dx = [,]x a b ∈ 0()y a y = 可以将区间[,]a b 分成n 段，那么方程在第i x 点有'()(,())i i i y x f x y x =，再用向前差商近似代替导数则为： ((1)()) (,())i i i i y x y x f x y x h +-= 在这里，h 是步长，即相邻两个结点间的距离。因此可以根据i x 点和i y 的数值计算出1i y +来： 1(,)i i i i y y h f x y +=+?0,1,2,i L = 这就是向前欧拉公式。 改进的欧拉公式：

将向前欧拉公式中的导数(,)i i f x y 改为微元两端导数的平均，即上式便是梯形的欧拉公式。 可见，上式是隐式格式，需要迭代求解。为了便于求解，使用改进的欧拉公式： 数值分析中，龙格－库塔法（Runge-Kutta ）是用于模拟常微分方程的解的重要的一类隐式或显式迭代法。实际上，龙格-库塔法是欧拉方法的一种推广，向前欧拉公式将导数项简单取为(,)n n f x y ,而改进的欧拉公式将导数项取为两端导数的平均。 龙格-库塔方法的基本思想： 在区间1[,]n n x x +内多取几个点，将他们的斜率加权平均，作为导数的近似。 龙格库塔法的家族中的一个成员如此常用，以至于经常被称为“RK4”或者就是“龙格库塔法”。 令初值问题表述如下。 '(,)y f t y =00()y t y = 则，对于该问题的RK4由如下方程给出： 11234(22)6 n n h y y k k k k +=++++ 其中 1(,)n n k f t y = 21(,)22 n n h h k f t y k =++ 32(,)22 n n h h k f t y k =++ 43(,)n n k f t h y hk =++

计算机仿真(龙格库塔方法)的软件VB设计与实现

计算机仿真（龙格库塔方法）的软件VB设计与实现 Dim a(0 To 10) As Single, b(0 To 10) As Single, c(0 To 10) As Single, d(0 To 10) As Single Dim e(0 To 10) As Single, h(0 To 10) As Single, p(0 To 10) As Single, q(0 To 10) As Single Dim n1 As Byte, n2 As Byte 'n1表示的是y 的阶数,n2表示的是输入函数的阶数 Dim i As Integer Dim f(0 To 2000) As Single 'f(m)=y,也就是各个时刻的y值 Dim X0(0 To 12) As Single, X1(0 To 12) As Single Dim dt As Single, u As Single 'dt为采样周期,u为输入Dim q1 As Single, p1 As Single, h1 As Single, e1 As Single Dim b0(0 To 10) As Single, bn(-4 To 10) As Single Dim m As Integer Dim max As Single Private Sub Combo1_Click() n1 = Combo1.ListIndex + 1 'n1表示的是y的阶数 For i = n1 + 1 To 8 Text1(i).Visible = False Label1(i).Visible = False Label9(i).Visible = False Next For i = 0 To n1 Text1(i).Visible = True Label1(i).Visible = True Label9(i).Visible = True Next Combo2.Clear Combo2.Text = "请选择输出u的最大阶数" For i = 0 To n1 Combo2.AddItem ((i) & "阶") Next End Sub Private Sub Combo2_Click() n2 = Combo2.ListIndex 'n2表示的是输入函数的阶数 For i = n2 + 1 To 8 Text2(i).Visible = False

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

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

可以得出四阶龙格－库塔公式，也就是在工程中应用广泛的经典龙格－库塔算法，即： 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); }

龙格库塔方法解题

龙格库塔方法 从常微分方程数值解法的几何意义看，欧拉方法取一点n t 处的斜率() n n Q t f k ,1 =作为平均斜率，因此欧拉方法近似公式为 t k Q Q n n ?+=+11 向后的欧拉方法则采用点1+n t 处的斜率 ()112,++=n n Q t f k 作为平均斜率，即 t k Q Q n n ?+=+21 所以这两种方法也称作矩形法。改进的欧拉方法则取点n t 处和点1+n t 处斜率 1k 和2k 的平均值作为平均斜率，即 ()t k k Q Q n n ?++=+2112 1 因此改进的欧拉方法又称为梯形方法。可以预见，若去多点处斜率的加权平均值作为平均斜率，误差会更小，这就是龙格-库塔方法。 最常用的是四阶龙格库卡近似计算公式，即： ()t k k k k Q Q n n ?++++=+43211226 1 式中 () () 314221312 121,2,2,,tk Q t f k k t Q t f k k t Q t f k Q t f k n n n n n n n n ?+=??? ? ???+=???? ???+==+++ 使用标准四阶龙格——库塔法求解初值问题

()(){}10210≤≤='=x xy y y 计算所用程序如下所示 #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++) { 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);

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 指导教师：日期： 教研室主任：日期：

四阶龙格库塔算法

//状态方程 /X=AX+BU化为微分方程求解 #include #include #include #include //定义运算步数； //定义步长； #define N 1000000 #define h 0.000001 double x[9]={1,2,3,4,5,6,7,8,9}; double u[4]={0}; //定义微分方程： double fx(double x[],double dx,int i) { //矩阵A和B double A[9][9]={{0,0,0,0,0,0,1,0.0023,0.0556},{0,0,0,0,0,0,0,0.9991,-0.0421},{0,0,0,0,0 ,0,0,0.0422,1.0007},{0.0001,-32.1478,0,-0.0197,0.004,0.0142,0.9141,1.7506,0.217} ,{32.1193,-0.0754,0,-0.0044,-0.0323,0.0046,-1.7475,0.9664,0.3464},{-1.3556,-1.78 91,0,0.0141,0.0031,-0.2634,0.1234,0.1345,-2.2448},{0,0,0,-0.0099,-0.031,0.0017,-3.2757,0.721,0.2452},{0,0,0,0.0044,-0.002,-0.0007,-0.1167,-0.6853,-0.0207},{0,0, 0,0.0009,0.0101,-0.0008,0.3469,-0.0784,-0.2875}}; double B[9][4]={{0,0,0,0},{0,0,0,0},{0,0,0,0},{0.0541,-0.0044,0.0218,-0.0007},{-0.0028, -0.0362,0.007,-0.0423},{0.0019,-0.0011,-0.3784,0},{0.001,-0.044,0.0154,-0.0339}, {-0.0136,0.001,-0.0013,-0.0001},{0,-0.004,-0.0162,0.0245}}; double a[9] = {0.0}; double b[9] = {0.0}; for(int j=0;j<9;j++) {a[i] += A[i][j]*(x[j]+dx);} for(int k=0;k<4;k++) {b[i] += B[k][0]*u[k];} return a[i]+b[i] ; } //主函数

各种龙格库塔计算

学生姓名：学号： 实验题目：数值积分法仿真应用 课程名称：控制系统仿真 实验目的： ?理解并掌握四种数值积分法的应用方法 ?理解并掌握应用数值积分法的Matlab程序设计 实验内容： ?应用MATLAB编写一个可以实现四种数值积分法运算的程序。 ?完善该程序使其可以用来对任意的传递函数进行仿真，可以实现其阶跃响 应和斜坡响应等，并能够对同一响应下不同积分法进行比较，要求与matlab 的Simulink仿真结果进行比较，同一积分法不同步长进行比较。 报告内容： （1）给出m文件的程序框图及其例题计算结果，并记录程序编写过程中出现的错误，并给出解决的方案。若没有得到解决，请说清楚你的问题。 （2）给出每一种情况的比较结果图形，并进行标注、说明。 此次实验约占整个科目成绩的25%，（其中程序部分占15%，报告部分占10%）提交日期：2011-4-22 实际提交日期： 声明：（包括报告内容和实验程序：此次实验完全自己完成；若得到了同学的帮助，请注明就哪一部分内容请教了同学，说明理由。） 签名： 注：本表打印出来作为封面，空格部分打印出来，手写填写姓名，学号和签名部分，实际提交日期由老师填写。 教师评语及成绩

题一：用4种方法求书本第86页第2-3题，输出响应y(t)在0≤t ≤1上，h=0.1时的数值解。1)0(,,=-=y y y 。要保留四位小数，并将结果与真解y(t)=e^(-t)比较，367879.01=-e 1. 欧拉法 程序： %欧拉法 T0=0; Tf=1; h=0.1; X=1;y=1;t=T0; N=round(Tf-T0)/h; for i=1:N f=-X; X=X+h*f; y=[y,X]; t=[t,t(i)+h]; end [t',y'] plot(t,y);grid on 运行结果如下：

欧拉法与龙格库塔法比较分析

解微分方程的欧拉法，龙格-库塔法简单实例比较 欧拉方法（Euler method)用以对给定初值的常微分方程(即初值问题)求解分为前EULER 法、后退EULER 法、改进的EULER 法。 缺点： 欧拉法简单地取切线的端点作为下一步的起点进行计算，当步数增多时，误差会因积累而越来越大。因此欧拉格式一般不用于实际计算。 改进欧拉格式（向前欧拉公式）： 为提高精度，需要在欧拉格式的基础上进行改进。采用区间两端的斜率的平均值作为直线方程的斜率。改进欧拉法的精度为二阶。 算法： 微分方程的本质特征是方程中含有导数项，数值解法的第一步就是设法消除其导数值。对于常微分方程： (,)dy f x y dx = [,]x a b ∈ 0()y a y = 可以将区间[,]a b 分成n 段，那么方程在第i x 点有'()(,())i i i y x f x y x =，再用向前差商近似代替导数则为： ((1)())(,())i i i i y x y x f x y x h +-= 在这里，h 是步长，即相邻两个结点间的距离。因此可以根据i x 点和i y 的数值计算出1i y +来： 1(,)i i i i y y h f x y +=+?0,1,2,i L = 这就是向前欧拉公式。 改进的欧拉公式：

将向前欧拉公式中的导数(,)i i f x y 改为微元两端导数的平均，即上式便是梯形的欧拉公式。 可见，上式是隐式格式，需要迭代求解。为了便于求解，使用改进的欧拉公式： 数值分析中，龙格－库塔法（Runge-Kutta ）是用于模拟常微分方程的解的重要的一类隐式或显式迭代法。实际上，龙格-库塔法是欧拉方法的一种推广，向前欧拉公式将导数项简单取为(,)n n f x y ,而改进的欧拉公式将导数项取为两端导数的平均。 龙格-库塔方法的基本思想： 在区间1[,]n n x x +内多取几个点，将他们的斜率加权平均，作为导数的近似。 龙格库塔法的家族中的一个成员如此常用，以至于经常被称为“RK4”或者就是“龙格库塔法”。 令初值问题表述如下。 '(,)y f t y =00()y t y = 则，对于该问题的RK4由如下方程给出： 11234(22)6 n n h y y k k k k +=++++ 其中 1(,)n n k f t y = 21(,)22 n n h h k f t y k =++ 32(,)22 n n h h k f t y k =++ 43(,)n n k f t h y hk =++

数值分析Matlab作业龙格库塔欧拉方法解二阶微分方程

数值分析M a t l a b作业龙格库塔欧拉方法解二阶微分方 程 -标准化文件发布号：（9456-EUATWK-MWUB-WUNN-INNUL-DDQTY-KII

Matlab 应用 使用Euler 和Rungkutta 方法解臂状摆的能量方程 背景 单摆是常见的物理模型，为了得到摆角θ的关于时间的函数，来描述单摆运动。由角动量定理我们知道 ε J M = 化简得到 0sin 22=+θθl g dt d 在一般的应用和计算中，只考虑摆角在5度以内的小摆动，因为可以吧简化为θ，这样比较容易解。实际上这是一个解二阶常微分方程的问题。 在这里的单摆是一种特别的单摆，具有均匀的质量M 分布在长为2的臂状摆上， 使用能量法建立方程 W T =

h mg ?=2J 2 1ω 化简得到 θθcos 35499.722=dt d 重力加速度取9.80665 1使用欧拉法 令dx dy z =，这样降阶就把二阶常微分方程转化为一阶微分方程组，再利用向前Euler 方法数值求解。 y(i+1)=y(i)+h*z(i); z(i+1)=z(i)+h*7.35499*cos(y(i)); y(0)=0 z(0)=0 精度随着h 的减小而更高，因为向前欧拉方法的整体截断误差与h 同阶，（因为是用了泰勒公式）所以欧拉方法的稳定区域并不大。

2.RK4-四阶龙格库塔方法 使用四级四阶经典显式Rungkutta公式 稳定性很好，RK4法是四阶方法，每步的误差是h5阶，而总积累误差为h4 阶。所以比欧拉稳定。 运行第三个程序：在一幅图中显示欧拉法和RK4法，随着截断误差的积累，欧 拉法产生了较大的误差 h=0.01