声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 4076|回复: 13

[混合编程] 关于用Runge-Kutta法求解非线性系统随机振动程序的问题

[复制链接]
发表于 2014-3-1 20:16 | 显示全部楼层 |阅读模式

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

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

x
本帖最后由 牛小贱 于 2014-3-4 15:30 编辑

我编了一段用Runge-Kutta法求解非线性系统随机振动的程序,不知道为什么老是出错,请各位大神帮忙看看!%%%%%%这是非线性系统的程序段
  1. function dx=yueh(t,x)
  2. ks=217582;m1=180;m2=2000;c1=4200;k2=2200000;c2=19540;
  3. dx=[x(2); (ks/m1)*x(3)-2*(c1/(2*m1))*2*(c1/(2*m1))*x(4)-(0.1*ks/m1)*(x(1)-x(3))^3-(ks/m1)*x(1);  x(4); c2*q1+k2*q+2*(c1/(2*m2))*x(2)-2*((c1+c2)/(2*m2))*x(4)+(ks/m2)*x(1)+(0.1*ks/m2)*(x(1)-x(3))^3-((k2+ks)/m2)*x(3)];
  4. %%%%%%%%%%主函数,suijisudu.txt和suijiweiyi.txt分别是随机振动的速度和位移
  5. clear,clc
  6. global q1 q;
  7. m=load('suijisudu.txt');
  8. n=load('suijiweiyi.txt');
  9. q1=m(:,2);
  10. q=n(:,2);
  11. [t,xE]=ode45(@yueh,[0,0.01,10],[0.1 0 0.1 0]);
  12. plot(t,xE(:,1));
复制代码




回复
分享到:

使用道具 举报

发表于 2014-3-2 16:56 | 显示全部楼层
请楼主上传一下suijisudu.txt和suijiweiyi.txt这两个文件!
 楼主| 发表于 2014-3-2 23:12 | 显示全部楼层
牛小贱 发表于 2014-3-2 16:56
请楼主上传一下suijisudu.txt和suijiweiyi.txt这两个文件!

不好意思,完了传了!现在补上!

suijisudu.txt

35.24 KB, 下载次数: 10

suijiweiyi.txt

36.38 KB, 下载次数: 9

发表于 2014-3-4 16:13 | 显示全部楼层
你可以参考这个帖子:http://www.ilovematlab.cn/thread-42808-1-1.html(关于Matlab中龙格-库塔(Runge-Kutta)方法原理及实现的 ),还有一个参考程序:调用函数:
  1. function xdot=odetest(t,x)
  2. %龙格-库塔算法测试
  3. %二阶非线性系统微分方程 xdotdot + 0.5*xdot + 2*x + x^2 = 0
  4. %求系统在初始条件为x(0)=1,xdot(0)=0的数值解
  5. xdot=zeros(2,1);
  6. xdot(1)=x(2);
  7. xdot(2)=-0.5*x(2)-2*x(1)-x(1)^2;
复制代码
主程序:


  1. tspan=[0 20];
  2. x0=[0 1];
  3. [t,x]=ode23(@odetest,tspan,x0);
  4. plot(t,x)
复制代码
对于你的这个问题,我觉得是这条语句有问题:dx=[x(2);(ks/m1)*x(3)-2*(c1/(2*m1))*2*(c1/(2*m1))*x(4)-(0.1*ks/m1)*(x(1)-x(3))^3-(ks/m1)*x(1);x(4);c2*q1+k2*q+2*(c1/(2*m2))*x(2)-2*((c1+c2)/(2*m2))*x(4)+(ks/m2)*x(1)+(0.1*ks/m2)*(x(1)-x(3))^3-((k2+ks)/m2)*x(3)];——
q1 和q都是1000*1的矩阵,显然运算有错。


点评

赞成: 5.0
赞成: 5
  发表于 2014-3-5 17:53
回复 支持 1 反对 0

使用道具 举报

发表于 2014-3-4 16:15 | 显示全部楼层
zhousidun 发表于 2014-3-2 23:12
不好意思,完了传了!现在补上!

你可以参考这个帖子:http://www.ilovematlab.cn/thread-42808-1-1.html(关于Matlab中龙格-库塔(Runge-Kutta)方法原理及实现的 ),还有一个参考程序:调用函数:
  1. function xdot=odetest(t,x)
  2. %龙格-库塔算法测试
  3. %二阶非线性系统微分方程 xdotdot + 0.5*xdot + 2*x + x^2 = 0
  4. %求系统在初始条件为x(0)=1,xdot(0)=0的数值解
  5. xdot=zeros(2,1);
  6. xdot(1)=x(2);
  7. xdot(2)=-0.5*x(2)-2*x(1)-x(1)^2;
复制代码
主程序:


  1. tspan=[0 20];
  2. x0=[0 1];
  3. [t,x]=ode23(@odetest,tspan,x0);
  4. plot(t,x)
复制代码
对于你的这个问题,我觉得是这条语句有问题:dx=[x(2);(ks/m1)*x(3)-2*(c1/(2*m1))*2*(c1/(2*m1))*x(4)-(0.1*ks/m1)*(x(1)-x(3))^3-(ks/m1)*x(1);x(4);c2*q1+k2*q+2*(c1/(2*m2))*x(2)-2*((c1+c2)/(2*m2))*x(4)+(ks/m2)*x(1)+(0.1*ks/m2)*(x(1)-x(3))^3-((k2+ks)/m2)*x(3)];——
q1 和q都是1000*1的矩阵,显然运算有错。


 楼主| 发表于 2014-3-4 16:50 | 显示全部楼层
牛小贱 发表于 2014-3-4 16:15
你可以参考这个帖子:http://www.ilovematlab.cn/thread-42808-1-1.html(关于Matlab中龙格-库塔(Runge-K ...

你的意思是应该把q和q1改成1X1000的矩阵么?
发表于 2014-3-4 17:19 | 显示全部楼层
zhousidun 发表于 2014-3-4 16:50
你的意思是应该把q和q1改成1X1000的矩阵么?

不是这个意思!因为c2*q1+k2*q+2*(c1/(2*m2))*x(2)-2*((c1+c2)/(2*m2))*x(4)+(ks/m2)*x(1)+(0.1*ks/m2)*(x(1)-x(3))^3-((k2+ks)/m2)*x(3) 包含着数与矩阵的运算,无法运算的,该式子前两个是个矩阵,后面全是赤裸裸的数,如何相加?我运行了你的程序,也改了一部分,但是就是有错误,一时半会找不出问题所在,所以才给你回帖子,给你写“资料”,希望你自己能从中找到错误之处,做出来也和我分享一下。你的程序,我再看看,改一下。

评分

1

查看全部评分

 楼主| 发表于 2014-3-4 20:42 | 显示全部楼层
牛小贱 发表于 2014-3-4 17:19
不是这个意思!因为c2*q1+k2*q+2*(c1/(2*m2))*x(2)-2*((c1+c2)/(2*m2))*x(4)+(ks/m2)*x(1)+(0.1*ks/m2)*( ...

感谢你细心的回答。其实我以前用过这个程序,只不过不是用随机数据的。那个时候我只是计算系统在不同周期信号下的幅频谱,程序是下面这个样子。虽然这个程序里面看起来是连续信号,但其实也是离散开了算的,我就想,如果我直接把数据离散开算,应该也能算。只是一直没成功。
function dx=yueh(t,x)
w=34.8;ks=217582;m1=180;m2=2000;c1=4200;k2=2200000;c2=19540;f=1500;
dx=[x(2); (ks/m1)*x(3)-2*(c1/(2*m1))*2*(c1/(2*m1))*x(4)-(0.1*ks/m1)*(x(1)-x(3))^3-(ks/m1)*x(1);  x(4); f*cos(w*t)+2*(c1/(2*m2))*x(2)-2*((c1+c2)/(2*m2))*x(4)+(ks/m2)*x(1)+(0.1*ks/m2)*(x(1)-x(3))^3-((k2+ks)/m2)*x(3)];
%%%%%%%%%%主函数,suijisudu.txt和suijiweiyi.txt分别是随机振动的速度和位移
clear,clc
[t,xE]=ode45(@yueh,[0,10],[0.1 0 0.1 0]);
plot(t,xE(:,1));
发表于 2014-12-29 11:52 | 显示全部楼层
学习了,我也想做一下随机振动。
发表于 2015-3-21 16:16 | 显示全部楼层
楼主确定非线性随机微分方程可以用MATLAB中的ode来求解吗,应该得先随机泰勒展开然后用龙哥库塔或欧拉等方法求解吧,如果楼主求解出来了,我想跟楼主讨教下方法,我也在解一个非线性随机微分方程组,犯愁好几天了
发表于 2015-3-26 20:35 | 显示全部楼层
ode求解常微分方程的吧。。。随机微分方程不可以吧?
发表于 2016-6-19 07:50 | 显示全部楼层
苏小贝 发表于 2015-3-26 20:35
ode求解常微分方程的吧。。。随机微分方程不可以吧?

感觉是啊,不能用ode吧,随机微分方程的话怎么求解呢?麻烦指教下啊
发表于 2016-6-20 08:58 | 显示全部楼层
  1. 试试。。。
  2. % o45ex1.m
  3. % solve y’’=tˆ2*y*cos(t)on [0,5] given y(0)=2, y’(0)=3
  4. % this can be written as a system of first order DEs
  5. % y1’=y2 y1(0)=2
  6. % y2’=tˆ2*y1*cos(t) y2(0)=3
  7. a=0; b=5; % interval
  8. y0=[2 3]; % initial condition
  9. tspan=[a:.1:b]; % interval of integration
  10. [t,y]=ode23(’o45ex1func’,tspan,y0);% soln in vectors t and y
  11. plot(t,y(:,1),t,y(:,2),’r--’) % plot y1 and y2 versus t
  12. legend(’y’,’dy/dt’,2)
  13. xlabel(’t’); ylabel(’y and dy/dt’)
  14. print -depsc o45ex1 % output to postscript file
  15. ———————————————————————————————————————
  16. function f=o45ex1func(t,y)
  17. % o45ex1func.m
  18. % the RHS of the system of DEs
  19. f(1)=y(2);
  20. f(2)=tˆ2.*y(1).*cos(t);
  21. f=f(:); % forces f to be a column vector
  22. return
复制代码


您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

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

GMT+8, 2024-11-28 08:33 , Processed in 0.082872 second(s), 23 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

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