随机微分方程的数值解法讨论！

无水1324

无聊，在这里想如下一个系统：
      d2x/dt2+e*dx/dt+x=f*sin(w*t)+f1                                           (1)
其中，e，f，w为常数，f1是一随机数，该随机数在0.7-1.3之内变化，且其均值为1。
没有学过随机微分系统，里面含有一个随机参数，暂且将方程(1)叫做随机微分方程，请高手指正。
问：
1、该类系统的分析可以用那些方法？
2、数值解法，能否像ode方程一样直接调用matlab函数求解呢？
3、matlab不行还有其它的数学计算工具行吗？

期待你的参与！

[ 本帖最后由 无水1324 于 2007-10-12 17:44 编辑 ]

wanyeqing2003

可以用线性加速度方法，非常简单。

无水1324

哦
好的谢谢！

请问直接用ode45会出现什么结果呢？

octopussheng

我是这么搞的，把随机数做成一个向量，然后用时变微分方程的ode求解方法计算，例子如下：
％ FUN2.M 带有数据集合的时变微分方程的例子
function  xdot = fun2 ( t, x, beta, omega, T, P )
pt = interp1 ( T, P, t );
xdot (2 ) = beta * x (2 ) – omega^2 * x (1) + pt ;
xdot (1 ) = x (2 );
xdot = xdot ( : );  % 使得xdot是一列
％ fun2.m 结束

在MATLAB中采用如下的形式调用这个函数：
beta = .1 ;
omega = 2 ;
A = .1 ;
w0 = 1.3 ;
theta = pi /4 ;
X0 = [ 0 1]’;
t0 = 0 ;
tf = 20 ;
T = t0 – eps : .1 : tf + theta +eps ;
P = A.* sin(w0 .*T –theta );
[ t, y ] = ode23 ( @ fun2 , [ t0, tf ], X0, [ ], theta, omega, T, P );
plot ( t, y (:,1));
hold on;
plot(t,y(:,2),'g');

无水1324

哦不错
pt = interp1 ( T, P, t );是什么意思？

octopussheng

这个就是用来产生时变项的命令！

无水1324

感谢，今天好好研究一下，:lol

无水1324

这个不对 ，我的是随机数，你这个是一些时变的数值.

octopussheng

其实随机数也可以直接生成的啊，然后变换到时变项就可以了

或者你可以用simulink生成时间－随机数的向量，然后用上面的方法同样可以实现的！

无水1324

我现在想，直接用ode45算会有什么问题？

octopussheng

直接算是不行的呀，还是需要作一些处理的！

无水1324

我就是在想为什么不可以，就是按照实际来考虑是可以的

octopussheng

数值计算时，随机数肯定是随时间变换的吧，这就形成一个时变的微分方程了，呵呵，又回到最开始的状态了吧！

无水1324

这个是对的，但是随机数我只知道其分布情况，所以你这个做法可能还是有点问题

octopussheng

通过随机数的分布即可以得到随机数随时间变化的向量，是上面的问题呀

无水觉得这种解决方法的问题出在哪里呢？赶紧讨论讨论，不然真错了我后面的东西不好做了！