MATLAB编程求教，谢谢回答

土木年华

各位大侠好， 麻烦问大家一个问题:

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';

土木年华

回复 1 # 土木年华 的帖子

请问大家inline这个内联函数的详细运用，当inline(fun), 函数fun如要有变化的参数怎么处理，谢谢大家，给指导一下！

ChaChing

先看看这帖, 看有没用
关于求解变参数微分方程 http://forum.vibunion.com/thread-44972-1-1.html

meiyongyuandeze

回复 1 # 土木年华 的帖子

想下你的是你的微分方程时什么样子的啊，怎么感觉既包含时间t又包含空间变量X，方便的话请上传个公式，看看，我帮你修改程序还是要看下公式！

雨人

一个小例子,希望对你有帮助！
tspan=[0 5];x0=5;
options=odeset('RelTol',1e-3);
% 内联函数求解
funs=inline('[2*x*u]','t','x','flag','u');
[t,y]=ode45(funs,tspan,x0,options,3);
plot(t,y)

saviour65

为什么在m文件中用sim调用.mdl文件没反应具体调用该怎么写代码 求大虾指导

土木年华

回复 4 # meiyongyuandeze 的帖子

你好， 首先感谢你的答复，我具体说一下我的问题，希望能得到你的帮助! 谢谢！！ 我这是个单自由度系统的振动方程：

方程描述

我就是想求出常微分方程的解(位移、速度、加速度），但是MATLAB水平不是很高，不大会处理具有变化参数的微分方程， 希望您给指导一下，谢谢

meiyongyuandeze

回复 7 # 土木年华 的帖子

我没有用龙格库塔，直接用的是ode45，你参考下吧！

1. function xdot=myfun(t,x)
   
2. mt=2+0.5*exp(-0.1*t)
   
3. kt=2+cos(t);
   
4. ct=1+0.22*sin(t);
   
5. xdot=[x(2);
   
6. -ct/mt*x(2)-kt/mt*x(1)];

复制代码

在命令窗口输入：

1. [t, y]=ode45('myfun',[0,100],[0,1])

复制代码

土木年华

回复 8 # meiyongyuandeze 的帖子

非常感谢你的指教， 我自己在做做，

土木年华

回复 8 # meiyongyuandeze 的帖子

不好意思，看了一下还是不大明白， 你这ode45输出的t和y代表什么，t是时间吧， y应该是速度吧，但如果我想同时输出位移和速度[x,y], 该怎么做呢， 不好意思，  麻烦你拉，谢谢你！！

meiyongyuandeze

回复 10 # 土木年华 的帖子

位移是y(:,1),速度y(:,2)

土木年华

回复 8 # meiyongyuandeze 的帖子

我想的是时间t变化为0:0.001:20s，初位移为10mm,初速度为10，求出后画出位移/时间曲线、速度/时间曲线 ， 非常感谢你！我弄了好长时间了一直没弄好出来，所以急着求你指教！

meiyongyuandeze

回复 12 # 土木年华 的帖子

1. [t, y]=ode45('myfun',[0,0.001:20],[0.01,10])
   

复制代码

土木年华

回复 13 # meiyongyuandeze 的帖子

3Q

土木年华

回复 8 # meiyongyuandeze 的帖子

你好！昨天得到你的知道很受启发，非常感谢你！
这几天我再网上一篇论文中学习了一个方法，一直在试着编程计算一下， 学习一下，想对我上面的这个例子运用这个方法计算一下，
我得出了位移时程曲线x(t),想利用下面方法求一个阻尼比

方法

最后就是想得出最后这个式子所描述的值，得到他的变化曲线，你能给我指导下吗？ matlab我也在学，但像这样的编程还是欠点火候，非常感谢！！