Search All of the Math Forum:
Views expressed in these public forums are not endorsed by
NCTM or The Math Forum.



Re: Stochastic differential equations with Simulink?
Posted:
Aug 29, 2013 8:15 AM


I'll try to answer myself. :)
I don't think Simulink can be used for SDDEs, and if(!) I got this right, the answer is so simple that Simulink wouldn't be the best choice anyway. In fact, I think this question is now more "how to solve SDDEs with Euler's method" than how to use Matlab to do this, but I'll post my code anyway to get some closure.. Also, I haven't verified the code below (I created it by combining bits and pieces of information on how to solve ODEs, SDEs, DDEs.. etc.), so if some SDDE wiz would take a look at it, it would be greatly appreciated.
So here is the code:
% clear all;close all;clc;
%%%%%%%%%%%%%%%%% %Parameters %%%%%%%%%%%%%%%%% TrialLength=100;%sec SampleFrequency=1000;%Hz dt=1/SampleFrequency;%sec TimeVec=(dt:dt:TrialLength); m=60;%kg g=9.81;%m/s2 h=1;%m J=60;%kgm2 Slope_a=0.4; P=0.25*m*g*h;%Nm/rad D=10;%Nms/rad K=0.8*m*g*h;%Nm/rad B=4;%Nms/rad TimeDelay=0.2;%s NoiseAmp=0.2;%Nm r=0.004;%radrad/s k=TimeDelay/dt;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Initial conditions and allocation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% X(1:TrialLength*SampleFrequency)=0;%Angle X(k+1)=0.01; V(1:TrialLength*SampleFrequency)=0;%Angle velocity
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Sway Angle with Euler's method: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for n=k+1:TrialLength*SampleFrequency+k a=((m*g*h)K)/J; b=B/J; c=P/J; d=D/J; e=NoiseAmp/J; N=randn(1,1);%Gaussian random number, mean=0, var=1 %Euler's method (hopefully correct..) %Angle X(n+1)=X(n)+V(n)*dt; %Angular velocity V(n+1)=V(n)+[a*X(n)+b*V(n)+c*X(nk)+d*V(nk)]*dt+e*N*sqrt(dt);
%PD active control on/off if X(nk)*(V(nk)Slope_a*X(nk))>0 && X(nk)^2+V(nk)^2>r^2 P=0.25*m*g*h; D=10; else P=0; D=0; end end
%Remove initial transient: X=X(k+2:end); V=V(k+2:end);
%Plotting figure;plot(X,V);ylabel('Angular velocity, (rad/sec)');xlabel('Angle (rad)'); figure;plot(TimeVec,X);ylabel('Angle (rad)');xlabel('Time (s)')
%
Aino



