rungkuta4.m

来自「4阶龙格--库塔法解常微分方程,比较好用,希望大家尝试.」· M 代码 · 共 21 行

function R=Rungkuta4(f,a,b,N,ya)
%  f是微分方程右端函数句柄
%  a,b自变量的取值范围[a,b]的左右端点。
%  N为区间等分个数或迭代步数。
%  h为步长。
%  ya为初值。
%  R=[x' y']是自变量x和y所组成的矩阵。
h=(b-a)/N;
x=zeros(1,N+1);
y=zeros(1,N+1);
x=a:h:b;
y(1)=ya;
for j=1:N
k1=feval(f,x(j),y(j));
k2=feval(f,x(j)+h/2,y(j)+h/2*k1);
k3=feval(f,x(j)+h/2,y(j)+h/2*k2);
k4=feval(f,x(j)+h,y(j)+h*k3);
y(j+1)=y(j)+(h/6)*(k1+2*k2+2*k3+k4);
end
R=[x' y']