方程
function dy=jeffcott(t,x)
dd=0.00016;
b=3;
% E=0.04;
cc=0.12;
f=0.12;
f0=25;
r=1.86;
e=sqrt(x(1)^2+x(3)^2);
dy=[x(2);
-2*cc*x(2)-x(1)+f*b/(1+b)*x(3)+b*dd*(x(1)-f*x(3))/[(1+b)*e]+r^2*cos(r*t);
x(4);
-2*cc*x(4)-x(3)-f*b/(1+b)*x(1)+b*dd*(f*x(1)+x(3))/[(1+b)*e]+r^2*sin(r*t)];
clc;
clear all;
x0=rand(4,1);
[t,y]=rk_41('jeffcott',[0,100,0.005],x0);
plot(t,y(:,1));
figure
plot(y(:,1),y(:,3)
rk_4程序
function [tout,yout]=rk_4(odefile,tspan,y0)
t0=tspan(1);
th=tspan(2);
if length(tspan)<=3,h=tspan(3);
else
h=tspan(2)-tsapn(1);
th=tspan(2);
end
tout=[t0:h:th]';
t=tout';
yout=[];
for t=tout'
k1=h*eval([odefile,'(t,y0)']);
k2=h*eval([odefile,'(t+h/2,y0+0.5*k1)']);
k3=h*eval([odefile,'(t+h/2,y0+0.5*k2)']);
k4=h*eval([odefile,'(t+h,y0+k3)']);
y0=y0+(k1+2*k2+2*k3+k4)/6;
yout=[yout;y0'];
end |