tqy0606 发表于 2010-11-3 16:25

求解微分方程组,方程组中带有零均值的白噪声序列(很急)

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

linlin820 发表于 2010-11-11 04:53

一个例子:
马尔可夫链近似求解随机微分方程程序的Matlab程序function mecca_zj
%马尔可夫链近似求解随机微分方程程序
%最简单的情况,利用随机模拟的思想直接模拟飞行航迹。
%利用多次模拟计算整体航迹空间分布的可能性

%构建空间网格%%%%%%%%%%%%%%%%%
M=500;
N=1000;
zMesh=repmat(0,);
Max=50;
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%模型参数和初值%%%%%%%%%%%%%%%%%%
sigma=1;
beita=1/5;
lamga=1/(4*(sigma^2));
X0=[-20;-30];
derta=0.1;
dertat=lamga*(derta^2);
timef=40;
stepMax=ceil(timef/dertat);
%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i=1:Max
%%%%%%%%%%%%%%%%%%%%%%%%%%%
step=1;
% 开始进入马尔可夫链的迭代过程
X=X0;
while(step<=stepMax)
    % X=pointP(:,step);
%   temp=(plan(step*dertat)+affine(,step*dertat,0)*X);%无风场
    temp=(plan(step*dertat)+affine(X,step*dertat,1));%存在仿射不变的风场
    ekq=temp(1)/(2*(sigma^2)*(1-hfunction(X,beita)));
    nkq=temp(2)/(2*(sigma^2)*(1-hfunction(X,beita)));
    xkq=1/(lamga*(sigma^2)*(1-hfunction(X,beita)))-4;
    ckq=2*cosh(-derta*ekq)+2*cosh(-derta*nkq)+xkq;
   
    pkl=exp(-derta*ekq)/ckq;
    pkr=exp(derta*ekq)/ckq;
    pkd=exp(-derta*nkq)/ckq;
    pku=exp(derta*nkq)/ckq;
    pko=xkq/ckq;
    %%%%%%%%%%%%%%%%%%%%%%%%%
   
    seed=rand(1);
    if(seed<pkl)
      X=X+[-derta;0];
    elseif(seed<pkl+pkr)
      X=X+;
    elseif(seed<pkl+pkr+pkd)
      X=X+;
    elseif(seed<pkl+pkr+pkd+pku)
      X=X+;
    else
      X=X;
    end
   
    J=round(N*(X(1)+30)/90);
    I=round(M*(X(2)+50)/80);
    if((I<=0)|(I>M))
      I=M;
    end
    if((J<=0)|(J>N))
      J=N;
    end
    zMesh(I,J)=zMesh(I,J)+1;
    step=step+1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%
i
end

temp=max(max(zMesh));
zMesh=255*zMesh/temp;
zMesh=uint8(zMesh);
imshow(zMesh);
imwrite(zMesh,'概率图.bmp','bmp');function v=plan(t)
%飞行器的飞行计划,且与位置无关
if(t<10)
    v=;
elseif(t<20)
    v=;
elseif(t<40)
    v=;
else
    v=;
endfunction f=affine(x,t,flag);
%风场的仿射变换矩阵函数

if(flag)
    R=(1/50)*;
    z=;
%   f=R*(x-z);
    f=R*x;
else
    f=0;
end
页: [1]
查看完整版本: 求解微分方程组,方程组中带有零均值的白噪声序列(很急)