补充一下，后面的那张图片就是那１８个公式的方程式，看起来可能会方便理解一点。
还有，我用runge-kutta四阶法求近似，
这里是runge-kutta的程序: rk4.m
function [x, y, dy] = rk4(deriv,n,x,dx,y)
x0 = x;
y0 = y;
[y,dy1] = feval(deriv,x0,y); % k1=f(x0,y0)
for i = 1:n
y(i) = y0(i) + 0.5*dx*dy1(i);
end
xm = x0 + 0.5*dx; % x0+0.5h
[y,dy2] = feval(deriv,xm,y);
% k2=f(x0+0.5h,y0+0.5k1*h)
for i = 1:n
y(i) = y0(i) + 0.5*dx*dy2(i); % x0+0.5*h, y0+0.5k2*h
end
[y,dy3] = feval(deriv,xm,y);
% k3=f(x0+0.5*h, y0+0.5k2*h)
for i = 1:n
y(i) = y0(i) + dx*dy3(i); % x0+h, y0+0.5k3*h
end
x = x0 + dx;
[y,dy] = feval(deriv,x,y);
% k4=f(x0+h,y0+0.5k3*h)
for i = 1:n
dy(i) = (dy1(i) + 2*(dy2(i) + dy3(i)) + dy(i))/6;
y(i) = y0(i) + dx*dy(i); % y
end
|