Search All of the Math Forum:

Views expressed in these public forums are not endorsed by NCTM or The Math Forum.

Notice: We are no longer accepting new posts, but the forums will continue to be readable.

Topic: Stochastic differential equations with Simulink?
Replies: 4   Last Post: Aug 29, 2013 8:15 AM

 Messages: [ Previous | Next ]
 Tikkuhirvi Tietavainen Posts: 99 Registered: 4/22/08
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;
TimeDelay=0.2;%s
NoiseAmp=0.2;%Nm
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(n-k)+d*V(n-k)]*dt+e*N*sqrt(dt);

%PD active control on/off
if X(n-k)*(V(n-k)-Slope_a*X(n-k))>0 && X(n-k)^2+V(n-k)^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

%-------------------------------------------------

-Aino

Date Subject Author
8/6/13 Tikkuhirvi Tietavainen
8/6/13 Phil Goddard
8/6/13 Phil Goddard
8/7/13 Tikkuhirvi Tietavainen
8/29/13 Tikkuhirvi Tietavainen