声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 1248|回复: 2

[非线性振动] 非线性系统的相平面的绘制

[复制链接]
发表于 2011-10-30 21:53 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

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




回复
分享到:

使用道具 举报

发表于 2011-10-31 14:17 | 显示全部楼层
除了有点乱 ,其他的还 好哦
 楼主| 发表于 2011-10-31 18:07 | 显示全部楼层
仅三个函数。。。不乱吧
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-5-22 17:18 , Processed in 0.121911 second(s), 17 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表