常微分方程組的四階RungeKutta龍格庫塔法matlab實(shí)現(xiàn)Word版

1、 常微分方程組的四階Runge-Kutta方法1.問題：1.1若用普通方法-僅適用于兩個方程組成的方程組編程實(shí)現(xiàn)：創(chuàng)建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為迭代次數(shù)%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

2、+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=

3、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)運(yùn)行結(jié)果：ans =1.0e+003 *0 3.0000 0.12000.0010 2.6637 0.

4、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數(shù)據(jù):情況二：對于x0=5000，y0=100命令行中輸入：>> rk4('f','g',0,10,5000,100,10)運(yùn)行結(jié)果：ans =1.0e+003 *0 5.0000 0.1

5、0000.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數(shù)據(jù)：結(jié)論：無論取得初值是哪一組，捕食者與被捕食者的數(shù)量總是一個增長另一個減少，并且是以T=5 為周期交替增長或減少的。這說明了捕食者與被捕食者是對立的，一個增長則另一個必定減少，從而達(dá)到食物鏈

6、相對穩(wěn)定。1.2改進(jìn)方法-一般方法代碼：創(chuàng)建M 文件：function x=rk4(f,x0,h,a,b,N)%x0初值%(a,b)為迭代區(qū)間%N迭代次數(shù)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以第一組數(shù)據(jù)為例：對于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)運(yùn)行結(jié)果：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 1