longgekuta.txt

％调用格式 =trapez_g( f_name ,a,b,n) ％f_name: 被积函数的文件名f(x) ％ a:x的上限 ％ b:x的下限 ％ 部分区间数 ％实例：trapez_g( s

a=0.4;c=0.2;gama=0.8;b=1.29;
h=0.01;
x(1)=0.1;y(1)=0;t(1)=0;n=0;
while t<=100
n=n+1;
k1=h*y(n);
m1=h*(-sin(x(n))-a*y(n)+b*cos(gama*t(n)).*sin(x(n))+c);
k2=h*(y(n)+m1/2);
m2=h*(-sin(x(n)+k1/2)-a*(y(n)+m1/2)+b*cos(gama*(t(n)+h/2)).*sin(x(n)+k1/2)+c);
k3=h*(y(n)+m2/2);
m3=h*(-sin(x(n)+k2/2)-a*(y(n)+m2/2)+b*cos(gama*(t(n)+h/2)).*sin(x(n)+k2/2)+c);
k4=h*(y(n)+m3);
m4=h*(-sin(x(n)+k3)-a*(y(n)+m3)+b*cos(gama*(t(n)+h)).*sin(x(n)+k3)+c);
x(n+1)=x(n)+(k1+2*k2+2*k3+k4)/6;
y(n+1)=y(n)+(m1+2*m2+2*m3+m4)/6;
t(n+1)=n*h;
end
figure(1)
plot(x(5000:end),y(5000:end),'r','MarkerSize',14)
title('phase of the warship equation')
figure(2)
plot(t(5000:end),y(5000:end),'r','MarkerSize',14)
title('Time of the warship equation')