声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 3918|回复: 20

[编程技巧] MATLAB编程求教,谢谢回答

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

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

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

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';
回复
分享到:

使用道具 举报

 楼主| 发表于 2011-4-9 20:16 | 显示全部楼层
回复 1 # 土木年华 的帖子

请问大家inline这个内联函数的详细运用,当inline(fun), 函数fun如要有变化的参数怎么处理,谢谢大家,给指导一下!
发表于 2011-4-10 00:58 | 显示全部楼层
先看看这帖, 看有没用
关于求解变参数微分方程 http://forum.vibunion.com/thread-44972-1-1.html

评分

1

查看全部评分

发表于 2011-4-10 11:47 | 显示全部楼层
回复 1 # 土木年华 的帖子

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

评分

1

查看全部评分

发表于 2011-4-11 08:18 | 显示全部楼层
一个小例子,希望对你有帮助!
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)

评分

1

查看全部评分

发表于 2011-4-12 11:53 | 显示全部楼层
为什么在m文件中用sim调用.mdl文件没反应具体调用该怎么写代码 求大虾指导
 楼主| 发表于 2011-4-14 16:06 | 显示全部楼层
本帖最后由 土木年华 于 2011-4-14 16:19 编辑

回复 4 # meiyongyuandeze 的帖子

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

方程描述

方程描述

我就是想求出常微分方程的解(位移、速度、加速度),但是MATLAB水平不是很高,不大会处理具有变化参数的微分方程, 希望您给指导一下,谢谢
发表于 2011-4-14 17:03 | 显示全部楼层
回复 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])
复制代码

评分

1

查看全部评分

 楼主| 发表于 2011-4-14 17:41 | 显示全部楼层
回复 8 # meiyongyuandeze 的帖子

非常感谢你的指教, 我自己在做做,
 楼主| 发表于 2011-4-14 18:04 | 显示全部楼层
回复 8 # meiyongyuandeze 的帖子

不好意思,看了一下还是不大明白, 你这ode45输出的t和y代表什么,t是时间吧, y应该是速度吧,但如果我想同时输出位移和速度[x,y], 该怎么做呢, 不好意思,  麻烦你拉,谢谢你!!
发表于 2011-4-14 18:13 | 显示全部楼层
回复 10 # 土木年华 的帖子

位移是y(:,1),速度y(:,2)
 楼主| 发表于 2011-4-14 18:20 | 显示全部楼层
回复 8 # meiyongyuandeze 的帖子

我想的是时间t变化为0:0.001:20s,初位移为10mm,初速度为10,求出后画出位移/时间曲线、速度/时间曲线 , 非常感谢你!我弄了好长时间了一直没弄好出来,所以急着求你指教!
发表于 2011-4-14 18:50 | 显示全部楼层
回复 12 # 土木年华 的帖子

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

评分

1

查看全部评分

 楼主| 发表于 2011-4-14 19:00 | 显示全部楼层
回复 13 # meiyongyuandeze 的帖子

3Q
 楼主| 发表于 2011-4-15 11:08 | 显示全部楼层
回复 8 # meiyongyuandeze 的帖子

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

方法

方法

最后就是想得出最后这个式子所描述的值,得到他的变化曲线,你能给我指导下吗? matlab我也在学,但像这样的编程还是欠点火候,非常感谢!!
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

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

GMT+8, 2024-11-28 17:20 , Processed in 0.065509 second(s), 23 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

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