求解微分方程组，方程组中带有零均值的白噪声序列（很急）

tqy0606

微分方程如下（状态方程）：
y1'=a*y1+b1*W(t)
y2'=c*y2+d*y1
其中W(t)为零均值的白噪声序列，求该方程在MATLAB中的解法，希望能给出详细的步骤，

linlin820

一个例子：
马尔可夫链近似求解随机微分方程程序的Matlab程序

1. function mecca_zj
   
2. %马尔可夫链近似求解随机微分方程程序
   
3. %最简单的情况,利用随机模拟的思想直接模拟飞行航迹。
   
4. %利用多次模拟计算整体航迹空间分布的可能性
   
5. 
   
6. %构建空间网格％％％％％％％％％％％％％％％％％
   
7. M=500;
   
8. N=1000;
   
9. zMesh=repmat(0,[M N]);
   
10. Max=50;
    
11. %%%%%%%%%%%%%%%%%%%%%%%%%%%
    
12. %模型参数和初值%%%%%%%%%%%%%%%%%%
    
13. sigma=1;
    
14. beita=1/5;
    
15. lamga=1/(4*(sigma^2));
    
16. X0=[-20;-30];
    
17. derta=0.1;
    
18. dertat=lamga*(derta^2);
    
19. timef=40;
    
20. stepMax=ceil(timef/dertat);
    
21. %%%%%%%%%%%%%%%%%%%%%%%%%%%
    
22. 
    
23. for i=1:Max
    
24. %%%%%%%%%%%%%%%%%%%%%%%%%%%
    
25. step=1;
    
26. % 开始进入马尔可夫链的迭代过程
    
27. X=X0;
    
28. while(step<=stepMax)
    
29.     % X=pointP(:,step);
    
30. %     temp=(plan(step*dertat)+affine([0;0],step*dertat,0)*X);％无风场
    
31.     temp=(plan(step*dertat)+affine(X,step*dertat,1));%存在仿射不变的风场
    
32.     ekq=temp(1)/(2*(sigma^2)*(1-hfunction(X,beita)));
    
33.     nkq=temp(2)/(2*(sigma^2)*(1-hfunction(X,beita)));
    
34.     xkq=1/(lamga*(sigma^2)*(1-hfunction(X,beita)))-4;
    
35.     ckq=2*cosh(-derta*ekq)+2*cosh(-derta*nkq)+xkq;
    
36.    
    
37.     pkl=exp(-derta*ekq)/ckq;
    
38.     pkr=exp(derta*ekq)/ckq;
    
39.     pkd=exp(-derta*nkq)/ckq;
    
40.     pku=exp(derta*nkq)/ckq;
    
41.     pko=xkq/ckq;
    
42.     %%%%%%%%%%%%%%%%%%%%%%%%%
    
43.    
    
44.     seed=rand(1);
    
45.     if(seed<pkl)
    
46.         X=X+[-derta;0];
    
47.     elseif(seed<pkl+pkr)
    
48.         X=X+[derta;0];
    
49.     elseif(seed<pkl+pkr+pkd)
    
50.         X=X+[0;-derta];
    
51.     elseif(seed<pkl+pkr+pkd+pku)
    
52.         X=X+[0;derta];
    
53.     else
    
54.         X=X;
    
55.     end
    
56.    
    
57.     J=round(N*(X(1)+30)/90);
    
58.     I=round(M*(X(2)+50)/80);
    
59.     if((I<=0)|(I>M))
    
60.         I=M;
    
61.     end
    
62.     if((J<=0)|(J>N))
    
63.         J=N;
    
64.     end
    
65.     zMesh(I,J)=zMesh(I,J)+1;
    
66.     step=step+1;
    
67. end
    
68. %%%%%%%%%%%%%%%%%%%%%%%%%%%
    
69. i
    
70. end
    
71. 
    
72. temp=max(max(zMesh));
    
73. zMesh=255*zMesh/temp;
    
74. zMesh=uint8(zMesh);
    
75. imshow(zMesh);
    
76. imwrite(zMesh,'概率图.bmp','bmp');

复制代码

1. function v=plan(t)
   
2. %飞行器的飞行计划,且与位置无关
   
3. if(t<10)
   
4.     v=[2;0];
   
5. elseif(t<20)
   
6.     v=[0;1];
   
7. elseif(t<40)
   
8.     v=[2;0];
   
9. else
   
10.     v=[0;0];
    
11. end

复制代码

1. function f=affine(x,t,flag);
   
2. %风场的仿射变换矩阵函数
   
3. 
   
4. if(flag)
   
5.     R=(1/50)*[0 1;-1 0];
   
6.     z=[3*t;(t.^2)/2];
   
7. %     f=R*(x-z);
   
8.     f=R*x;
   
9. else
   
10.     f=0;
    
11. end

复制代码