资源预览内容
第1页 / 共5页
第2页 / 共5页
第3页 / 共5页
第4页 / 共5页
第5页 / 共5页
亲,该文档总共5页全部预览完了,如果喜欢就下载吧!
资源描述
常微分方程组的四阶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:Nf1=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.12000.0010 2.6637 0.09260.0020 3.7120 0.07740.0030 5.5033 0.08860.0040 4.9866 0.11930.0050 3.1930 0.11950.0060 2.7665 0.09510.0070 3.6543 0.07990.0080 5.2582 0.08840.0090 4.9942 0.11570.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.10000.0010 4.1883 0.11440.0020 3.2978 0.10720.0030 3.3468 0.09220.0040 4.2020 0.08760.0050 4.8807 0.09950.0060 4.2090 0.11260.0070 3.3874 0.10690.0080 3.4011 0.09340.0090 4.1568 0.08890.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+1L1=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 93.0000 2.6637 3.7120 5.5033 4.9866 3.1930 2.7665 3.6543 5.25820.1200 0.0926 0.0774 0.0886 0.1193 0.1195 0.0951 0.0799 0.0884Columns 10 through 124.9942 3.3541 2.86890.1157 0.1185 0.0970结果分析:与1.1 的结果相同,当然这个方法还适合于一般方程组。
网站客服QQ:2055934822
金锄头文库版权所有
经营许可证:蜀ICP备13022795号 | 川公网安备 51140202000112号