2楼
function [ x,y] = nark4v( dyfun,xspan,y0,e,h )
%用途：变步长4阶经典Runge-Kutta格式解常微分方程y'=f(x,y),y(xo)=y0
%格式：[x,y]=nark4(defun,xspan,y0,h) dyfun为函数f(x,y),xspan为求解区间
% [x0,xN],y0为初值y（x0），h为初始步长（默认为xspan的1/10），x返回节点，
% y返回数值解，e为精度要求
%feval为函数求值命令，'为转置，nargin为输入变量个数,eps为浮点数识别精度2^（-52）
if nargin<5,h=(xspan(2)-xspan(1))/10;end
n=1;x(n)=xspan(1);y(n)=y0;
[y1,y2]=comput(dyfun,x(n),y(n),h);
while x(n)
if abs(y2-y1)/10>e
while abs(y2-y1)/10>e
h=h/2;
[y1,y2]=comput(dyfun,x(n),y(n),h);
end
else
while abs(y2-y1)/10<=e
h=2*h;
[y1,y2]=comput(dyfun,x(n),y(n),h);
end
h=h/2;h=min(h,xspan(2)-x(n));
[y1,y2]=comput(dyfun,x(n),y(n),h);
end
n=n+1;x(n)=x(n-1)+h;y(n)=y2;
[y1,y2]=comput(dyfun,x(n),y(n),h);
end
x=x';y=y';
function [y1,y2]=comput(dyfun,x,y,h)
y1=rk4(dyfun,x,y,h);
y21=rk4(dyfun,x,y,h/2);
y2=rk4(dyfun,x+h/2,y21,h/2);
function y=rk4(dyfun,x,y,h)
k1=feval(dyfun,x,y);
k2=feval(dyfun,x+h/2,y+h/2*k1);
k3=feval(dyfun,x+h/2,y+h/2*k2);
k4=feval(dyfun,x+h,y+h*k3);
y=y+h*(k1+2*k2+2*k3+k4)/6;
回复
