文档库 最新最全的文档下载
当前位置:文档库 › 常微分方程组的四阶RungeKutta龙格库塔法matlab实现

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

常微分方程组的四阶RungeKutta龙格库塔法matlab实现
常微分方程组的四阶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

数据:

情况二:对于x0=5000,y0=100

命令行中输入:

>> rk4('f','g',0,10,5000,100,10)

运行结果:

ans =

1.0e+003 *

0 5.0000 0.1000

0.0010 4.1883 0.1144

0.0020 3.2978 0.1072

0.0030 3.3468 0.0922

0.0040 4.2020 0.0876

0.0050 4.8807 0.0995

0.0060 4.2090 0.1126

0.0070 3.3874 0.1069

0.0080 3.4011 0.0934

0.0090 4.1568 0.0889

0.0100 4.7753 0.0991

数据:

结论:无论取得初值是哪一组,捕食者与被捕食者的数量总是一个增长另一个减少,并且是以T=5 为周期交替增长或减少的。这说明了捕食者与被捕食者是对立的,一个增长则另一个必定减少,从而达到食物链相对稳定。

1.2改进方法---------一般方法

代码:

创建M 文件:

function x=rk4(f,x0,h,a,b,N)

%x0初值

%(a,b)为迭代区间

%N迭代次数

t=a:h:b;

t=zeros(1:N+1);

x(:,1)=x0;

for i=1:N+1

L1=f(t(i),x(:,i));

L2=f(t(i)+h/2,x(:,i)'+(h/2)*L1);

L3=f(t(i)+h/2,x(:,i)'+(h/2)*L2);

L4=f(t(i)+h,x(:,i)'+h*L3);

x(:,i+1)=x(:,i)'+(h/6)*(L1+2*L2+2*L3+L4);

end

以第一组数据为例:

对于x0=3000,y0=120

控制台中输入:

f=@(t,x)[2*x(1)-0.02*x(1)*x(2),0.0002*x(1)*x(2)-0.8*x(2)]; x0=[3000,120];

N=10;

a=0;

b=10;

rk4(f,x0,a,b,10)

运行结果:

ans =

1.0e+003 *

Columns 1 through 9

3.0000 2.6637 3.7120 5.5033

4.9866 3.1930 2.7665 3.6543

5.2582 0.1200 0.0926 0.0774 0.0886 0.1193 0.1195 0.0951 0.0799 0.0884 Columns 10 through 12

4.9942 3.3541 2.8689

0.1157 0.1185 0.0970

结果分析:与1.1 的结果相同,当然这个方法还适合于一般方程组。

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实现~ 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 四阶龙格-库塔法求微分方程

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文本文件)运行;

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

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.实验目的 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

龙格库塔法求微分方程2

《MATLAB 程序设计实践》课程考核 一、编程实现“四阶龙格-库塔(R-K )方法求常微分方程”,并举一 例应用之。 【实例】采用龙格-库塔法求微分方程: ?? ?==+-=0 , 0)(1 '00 x x y y y 1、算法说明: 在龙格-库塔法中,四阶龙格-库塔法的局部截断误差约为o(h5),被广泛应用于解微分方程的初值问题。其算法公式为: )22(6 3211k k k h y y n n +++=+ 其中: ?????????++=++=++ ==) ,() 21 ,21()21 ,21() ,(34 23121hk y h x f k hk y h x f k hk y h x f k y x f k n n n n n n n n 2、流程图: 2.1、四阶龙格-库塔(R-K )方法流程图:

2.2、实例求解流程图:

3、源程序代码 3.1、四阶龙格-库塔(R-K)方法源程序: function [x,y] = MyRunge_Kutta(fun,x0,xt,y0,PointNum,varargin) %Runge-Kutta 方法解微分方程形为 y'(t)=f(x,y(x)) %此程序可解高阶的微分方程。只要将其形式写为上述微分方程的向量形式 %函数 f(x,y): fun %自变量的初值和终值:x0, xt %y0表示函数在x0处的值,输入初值为列向量形式 %自变量在[x0,xt]上取的点数:PointNum %varargin为可输入项,可传适当参数给函数f(x,y) %x:所取的点的x值 %y:对应点上的函数值 if nargin<4 | PointNum<=0 PointNum=100; end if nargin<3 y0=0; end y(1,:)=y0(:)'; %初值存为行向量形式h=(xt-x0)/(PointNum-1); %计算步长 x=x0+[0:(PointNum-1)]'*h; %得x向量值 for k=1:(PointNum)%迭代计算 f1=h*feval(fun,x(k),y(k,:),varargin{:}); f1=f1(:)'; %得公式k1 f2=h*feval(fun,x(k)+h/2,y(k,:)+f1/2,varargin{:}); f2=f2(:)'; %得公式k2 f3=h*feval(fun,x(k)+h/2,y(k,:)+f2/2,varargin{:}); f3=f3(:)'; %得公式k3 f4=h*feval(fun,x(k)+h,y(k,:)+f3,varargin{:}); f4=f4(:)'; %得公式k4 y(k+1,:)=y(k,:)+(f1+2*(f2+f3)+f4)/6; %得y(n+1) end 3.2、实例求解源程序: %运行四阶R-K法

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

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

例题一 四阶龙格-库塔法的具体形式为: 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

基于matlab的龙格库塔法求解布拉修斯方程

Runge —Kutta 法求解布拉修斯解 摘要 薄剪切层方程主要有三种解法,即相似解,非相似条件下对偏微分方程组的数值解和近似解。布拉修斯解是布拉修斯于1908年求出的,它是零攻角沿平板流动的相似解。本文用四阶Runge —Kutta 法求解高阶微分方程的方法,并用matlab 编程实现,求得了与实际层流边界层相符合的数值解。 关键词:布拉修斯解,相似解,Runge —Kutta 法,数值解。 1 布拉修斯近似解方程 二维定常不可压缩层流边界层的方程为: =??+??y v x u (1) 2 2 y u v dx d y u v x u u u u e e ??+= ??+?? (2) 边界条件为 :0=y )(,0x v u v w = = :δ=y )(x u u e = 将式(1)和式(2)进行法沃克纳—斯坎变换(简称F —S 变换),将边界层方程无量纲化,即设 y x v u e 5 .0?? ? ? ??=η (3) x x = (4) 得出F —S 变换后的动量方程 () []()[] ??? ? ???''-?'?'='-+''++'''+x f f x f f x f m f f m f t k 2 21211 (5) 其中k 为流动类型指标,横曲率项t 为

2 1 2 1 2 0c o s 211???? ????? ? ??? ? ????? ? ??+ +-=ηφ e u vx L r L t (6) m 是量纲一的压力梯度参数,定义为 x d du u x m e e = (7) 其边界条件变为 :0=η 0='f :∞=η 1='f 对于二维平面实壁流动(:0=η0=w f )可以忽略横曲率项t 的轴对称流动,式(5)成为 () [ ]?? ? ? ???' '-?'?'='-+''++ '''x f f x f f x f m f f m f 2 12 1 (8) 根据相似解的定义,方程(8)中的函数f 若式相似的,则它应只与η有关而与x 无关,即对x 的偏导数应为零。于是方程(8)应成为 () [ ]012 12 ='-+''++ '''f m f f m f (9) 若f w 为常数,则方程(9)的边界条件为 :0=η 常数==w f f ; 0='='w f f :∞=η 1='f 2 布拉修斯解 布拉修斯于1908年求出了零攻角沿平板流动的解。这时 0==m u e 常数; 因而方程(9)成为 02 1=''+ '''f f f (10) 此即布拉修斯方程。对于实壁,0=w f ,边界条件成为 :0=η 0==w f f ; 0='='w f f :∞=η 1='f

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

四阶龙格库塔法原理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

标准四阶龙格——库塔法

实验名:常微分方程数值解法 实习目的: (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 展开式有较多的项达到一致,从而得出较高阶的数值公式,这就是龙格—库塔法的基本思想。

常微分方程组的四阶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 数据:

龙格-库塔法(Runge-Kutta)matlab代码及含义

龙格-库塔法(Runge-Kutta) 数值分析中,龙格-库塔法(Runge-Kutta)是用于模拟常微分方程的解的重要的一类隐式或显式迭代法。这些技术由数学家卡尔·龙格和马丁·威尔海姆·库塔于1900年左右发明。 经典四阶龙格库塔法 RK4””或者就是龙格库塔法的家族中的一个成员如此常用,以至于经常被称为“RK4 “龙格库塔法”。 令初值问题表述如下。 则,对于该问题的RK4由如下方程给出: 其中 这样,下一个值(yn+1)由现在的值(yn)加上时间间隔(h)和一个估算的斜率的乘积决定。该斜率是以下斜率的加权平均: k1是时间段开始时的斜率; k2是时间段中点的斜率,通过欧拉法采用斜率k1来决定y在点tn+h/2的值;k3也是中点的斜率,但是这次采用斜率k2决定y值; k4是时间段终点的斜率,其y值用k3决定。 当四个斜率取平均时,中点的斜率有更大的权值: RK4法是四阶方法,也就是说每步的误差是h5阶,而总积累误差为h4阶。 注意上述公式对于标量或者向量函数(y可以是向量)都适用。 显式龙格库塔法 显示龙格-库塔法是上述RK4法的一个推广。它由下式给出 其中

(注意:上述方程在不同著述中由不同但却等价的定义)。 要给定一个特定的方法,必须提供整数s(阶段数),以及系数aij(对于1≤j

龙格库塔法RKF45Matlab实现

龙格库塔法RKF45的Matlab实现 2007-08-16 14:03:32| 分类:MatLab/Maple/Mat|字号订阅 4阶5级龙格库塔法用于解一阶微分方程(组),对于高阶微分方程,可以将其转换为一阶微分方程组求解。原程序由John.H.Mathews编写(数值方法matlab版),但只能解微分方程,不能解微分方程组。由LiuLiu@uestc修改,使之能够解微分方程组。该程序精度比matlab自带的ode45更高。 rkf45.m: function [Rt Rx]=rkf45(f,tspan,ya,m,tol) % Input: % - f function column vector % - tspan[a,b] left & right point of [a,b] % - ya initial value column vector % -m initial guess for number of steps % -tol tolerance % Output: % - Rt solution: vector of abscissas % - Rx solution: vector of ordinates % Program by John.Mathews, improved by liuliu@uestc if length(tspan)~=2 error('length of vector tspan must be 2.'); end if ~isnumeric(tspan) error('TSPAN should be a vector of integration steps.'); end if ~isnumeric(ya) error('Ya should be a vector of initial conditions.'); end h = diff(tspan); if any(sign(h(1))*h <= 0) error('Entries of TSPAN are not in order.') ; end a=tspan(1);

单步龙格库塔比例导引弹道计算matlab源程序

单步龙格库塔比例导引弹道计算matlab源程序 本程序特别适合于弹道计算等方面使用。比ODE45函数的速度加快了很多,且程序充分展示了Matlab向量运算的强大功能,以及编程的简单快捷。 但愿本程序,可以为大家提高工作效率提高帮助。 一、rightF.m文件(调用的时候,要根据自己的应用,改变右函数值) %右函数定义的,用于单步龙格库塔法计算。 %如果需要使用全局变量,要使用global定义。 %这种格式的,也完全可以被ode45函数调用。 function dy=rightF(t,y)%right functions dy=zeros(2,1);%a column vector dy(1)=y(2); dy(2)=3*y(2)/t-4*y(1)/(t^2); 二、stepRK.m文件(本函数任何情况都不用变的) %函数说明,t表示自变量(一般为时间),y表示变量的初值(是一个行向量),h表示步长 %运行的结果,也是一个行向量 function dy=stepRK(t,y,h) K1=rightF(t,y); K2=rightF(t+h/2,y'+h*K1/2); K3=rightF(t+h/2,y'+h*K2/2); K4=rightF(t+h,y'+h*K3); dy=y'+(K1+2*K2+2*K3+K4)*h/6; dy=dy'; 三、示例程序(rk_test.m) %龙格库塔法测试程序,与欧阳联渊书一致 t=1.; a=[10]; h=0.1; for i=1:10 b=stepRK(t,a,h) t=t+h; a=b; arrayT(i)=i; array1(i)=b(1); end plot(arrayT,array1,'*')

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

相关文档
相关文档 最新文档