求解微分方程组，方程组中带有零均值的白噪声序列（很急）
微分方程如下（状态方程）：y1'=a*y1+b1*W(t)
y2'=c*y2+d*y1
其中W(t)为零均值的白噪声序列，求该方程在MATLAB中的解法，希望能给出详细的步骤，
一个例子：
马尔可夫链近似求解随机微分方程程序的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]