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

zhousidun

我编了一段用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));

复制代码

牛小贱

请楼主上传一下suijisudu.txt和suijiweiyi.txt这两个文件！

zhousidun

牛小贱 发表于 2014-3-2 16:56
请楼主上传一下suijisudu.txt和suijiweiyi.txt这两个文件！

不好意思，完了传了！现在补上！

牛小贱

你可以参考这个帖子：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的矩阵，显然运算有错。

牛小贱

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的矩阵，显然运算有错。

zhousidun

牛小贱 发表于 2014-3-4 16:15
你可以参考这个帖子：http://www.ilovematlab.cn/thread-42808-1-1.html（关于Matlab中龙格-库塔(Runge-K ...

你的意思是应该把q和q1改成1X1000的矩阵么？

牛小贱

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) 包含着数与矩阵的运算，无法运算的，该式子前两个是个矩阵，后面全是赤裸裸的数，如何相加？我运行了你的程序，也改了一部分，但是就是有错误，一时半会找不出问题所在，所以才给你回帖子，给你写“资料”，希望你自己能从中找到错误之处，做出来也和我分享一下。你的程序，我再看看，改一下。

zhousidun

牛小贱 发表于 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));

dollfish000

学习了，我也想做一下随机振动。

苏小贝

楼主确定非线性随机微分方程可以用MATLAB中的ode来求解吗，应该得先随机泰勒展开然后用龙哥库塔或欧拉等方法求解吧，如果楼主求解出来了，我想跟楼主讨教下方法，我也在解一个非线性随机微分方程组，犯愁好几天了

苏小贝

ode求解常微分方程的吧。。。随机微分方程不可以吧？

lucky211

苏小贝 发表于 2015-3-26 20:35
ode求解常微分方程的吧。。。随机微分方程不可以吧？

感觉是啊，不能用ode吧，随机微分方程的话怎么求解呢？麻烦指教下啊

敷衍会致命

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
    

复制代码