马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
%要求的一阶方程组的形式
%dx=f(t,x,y); %位移的导数,x代表位移,y代表速度
%dy=g(t,x,y); %速度的导数
%在末尾出设置
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%二阶常微分方程组的求解——画相位图
function Runge_Kutta()
clc
clear all
%算法参数预设
h=0.1;Tmax=0.6; %步长为h,计算时间为20s,n用于计数
xold=0;yold=2; %设置初始值
%四阶龙格算法
for t=0:h:Tmax
t=0;
m1=f(t,xold,yold)*h;k1=g(t,xold,yold)*h; %将h放在这里减少计算舍去误差
m2=f(t+0.5*h,xold+0.5*m1,yold+0.5*k1)*h;k2=g(t+0.5*h,xold+0.5*m1,yold+0.5*k1)*h;
m3=f(t+0.5*h,xold+0.5*m2,yold+0.5*k2)*h;k3=g(t+0.5*h,xold+0.5*m2,yold+0.5*k2)*h;
m4=f(t+h,xold+m3,yold+k3)*h;k4=g(t+h,xold+m3,yold+k3)*h;
%保存
if t==0; n=1;end
TT(n,:)=t;XX(n,:)=xold;YY(n,:)=yold;n=n+1;
%更新
xold=xold+(m1+2*m2+2*m3+m4)/6;
yold=yold+(k1+2*k2+2*k3+k4)/6;
t=t+h;
end
%计算结果输出
%plot(TT,XX,'b');grid on;
plot(XX,YY,'r');grid on
end
%要求的位移的导数表达式dx=f(t,x,y);
function dx=f(t,x,y)
dx=2*x+4*y; %不同的系统,需要修改
end
%要求的速度的导数表达式dy=g(t,x,y);
function dy=g(t,x,y)
%t是时间,x是位移,或者速度的导数,y是速度,或者速度的导数
dy=-x+6*y; %不同的系统,需要修改
end
|