关于用Runge-Kutta法求解非线性系统随机振动程序的问题
本帖最后由 牛小贱 于 2014-3-4 15:30 编辑我编了一段用Runge-Kutta法求解非线性系统随机振动的程序，不知道为什么老是出错，请各位大神帮忙看看！%%%%%%这是非线性系统的程序段
function dx=yueh(t,x)
ks=217582;m1=180;m2=2000;c1=4200;k2=2200000;c2=19540;
dx=;
%%%%%%%%%%主函数，suijisudu.txt和suijiweiyi.txt分别是随机振动的速度和位移
clear,clc
global q1 q;
m=load('suijisudu.txt');
n=load('suijiweiyi.txt');
q1=m(:,2);
q=n(:,2);
=ode45(@yueh,,);
plot(t,xE(:,1));
请楼主上传一下suijisudu.txt和suijiweiyi.txt这两个文件！ 牛小贱 发表于 2014-3-2 16:56
请楼主上传一下suijisudu.txt和suijiweiyi.txt这两个文件！
不好意思，完了传了！现在补上！
你可以参考这个帖子：http://www.ilovematlab.cn/thread-42808-1-1.html（关于Matlab中龙格-库塔(Runge-Kutta)方法原理及实现的 ），还有一个参考程序：调用函数：function xdot=odetest(t,x)
%龙格-库塔算法测试
%二阶非线性系统微分方程 xdotdot + 0.5*xdot + 2*x + x^2 = 0
%求系统在初始条件为x(0)=1,xdot(0)=0的数值解
xdot=zeros(2,1);
xdot(1)=x(2);
xdot(2)=-0.5*x(2)-2*x(1)-x(1)^2;
主程序：
tspan=;
x0=;
=ode23(@odetest,tspan,x0);
plot(t,x)对于你的这个问题，我觉得是这条语句有问题：dx=;——q1 和q都是1000*1的矩阵，显然运算有错。
zhousidun 发表于 2014-3-2 23:12
不好意思，完了传了！现在补上！
你可以参考这个帖子：http://www.ilovematlab.cn/thread-42808-1-1.html（关于Matlab中龙格-库塔(Runge-Kutta)方法原理及实现的 ），还有一个参考程序：调用函数：function xdot=odetest(t,x)
%龙格-库塔算法测试
%二阶非线性系统微分方程 xdotdot + 0.5*xdot + 2*x + x^2 = 0
%求系统在初始条件为x(0)=1,xdot(0)=0的数值解
xdot=zeros(2,1);
xdot(1)=x(2);
xdot(2)=-0.5*x(2)-2*x(1)-x(1)^2;
主程序：
tspan=;
x0=;
=ode23(@odetest,tspan,x0);
plot(t,x)对于你的这个问题，我觉得是这条语句有问题：dx=;——q1 和q都是1000*1的矩阵，显然运算有错。
牛小贱 发表于 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) 包含着数与矩阵的运算，无法运算的，该式子前两个是个矩阵，后面全是赤裸裸的数，如何相加？我运行了你的程序，也改了一部分，但是就是有错误，一时半会找不出问题所在，所以才给你回帖子，给你写“资料”，希望你自己能从中找到错误之处，做出来也和我分享一下{:{04}:}。你的程序，我再看看，改一下。 牛小贱 发表于 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=;
%%%%%%%%%%主函数，suijisudu.txt和suijiweiyi.txt分别是随机振动的速度和位移
clear,clc
=ode45(@yueh,,);
plot(t,xE(:,1)); 学习了，我也想做一下随机振动。 楼主确定非线性随机微分方程可以用MATLAB中的ode来求解吗，应该得先随机泰勒展开然后用龙哥库塔或欧拉等方法求解吧，如果楼主求解出来了，我想跟楼主讨教下方法，我也在解一个非线性随机微分方程组，犯愁好几天了 ode求解常微分方程的吧。。。随机微分方程不可以吧？ 苏小贝 发表于 2015-3-26 20:35
ode求解常微分方程的吧。。。随机微分方程不可以吧？
感觉是啊，不能用ode吧，随机微分方程的话怎么求解呢？麻烦指教下啊 试试。。。
% o45ex1.m
% solve y’’=tˆ2*y*cos(t)on given y(0)=2, y’(0)=3
% this can be written as a system of first order DEs
% y1’=y2 y1(0)=2
% y2’=tˆ2*y1*cos(t) y2(0)=3
a=0; b=5; % interval
y0=; % initial condition
tspan=; % interval of integration
=ode23(’o45ex1func’,tspan,y0);% soln in vectors t and y
plot(t,y(:,1),t,y(:,2),’r--’) % plot y1 and y2 versus t
legend(’y’,’dy/dt’,2)
xlabel(’t’); ylabel(’y and dy/dt’)
print -depsc o45ex1 % output to postscript file
———————————————————————————————————————
function f=o45ex1func(t,y)
% o45ex1func.m
% the RHS of the system of DEs
f(1)=y(2);
f(2)=tˆ2.*y(1).*cos(t);
f=f(:); % forces f to be a column vector
return
页:
[1]