马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
各位大侠好, 麻烦问大家一个问题:
t=0:0.001:20;
mt=(2+exp(-0.1*t));
kt=(2+cos(t));
ct=(1+sin(t));
wt=zeros(20001,1);
qt=zeros(20001,1);
for k=1:20001;
wt(k)=(1/(2*pi))*sqrt(kt(k)/mt(k));
qt(k)=ct(k)/(2*sqrt(mt(k)*kt(k)));
dyfun=inline('(ct(k)/mt(k))*y+(kt(k)/mt(k))*x'); %y为速度,x为位移
[x,y]=nark4(dyfun,[-10,10],10,0.001); %由四阶龙格库塔公式求解
k=k+1;
end
我这个就是想用四阶龙格库塔公式求解微分方程, 但微分方程里边的参数mt,ct,kt是变化的, 但我用个循环去做,但怎么运行不出来,哪位大侠做过,给指导一下,谢谢!!
function[x,y]=nark4(dyfun,xspan,y0,h)
% 4阶经典Runge-Kutta格式解常微分方程y′=f(x,y),y(x0)=y0
% dyfun :函数f(x,y)
% xspan: 求解区间[x0,xN]
% h;步长
% y0: 初值y(x0)
% x返回节点, y 返回数值解
x=xspan(1):h:xspan(2);
y(1)=y0;
for n=1:length(x)-1
k1=feval(dyfun,x(n),y(n));
k2=feval(dyfun,x(n)+h/2,y(n)+h/2*k1);
k3=feval(dyfun,x(n)+h/2,y(n)+h/2*k2);
k4=feval(dyfun,x(n+1),y(n)+h*k3);
y(n+1)=y(n)+h*(k1+2*k2+2*k3+k4)/6;
end
x=x';y=y';
|