Drexel dragonThe Math ForumDonate to the Math Forum



Search All of the Math Forum:

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


Math Forum » Discussions » Software » comp.soft-sys.matlab

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

Advanced Search

Back to Topic List Back to Topic List Jump to Tree View Jump to Tree View   Messages: [ Previous | Next ]
Tikkuhirvi Tietavainen

Posts: 99
Registered: 4/22/08
Re: Stochastic differential equations with Simulink?
Posted: Aug 29, 2013 8:15 AM
  Click to see the message monospaced in plain text Plain Text   Click to reply to this topic Reply

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;%rad-rad/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(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
figure;plot(X,V);ylabel('Angular velocity, (rad/sec)');xlabel('Angle (rad)');
figure;plot(TimeVec,X);ylabel('Angle (rad)');xlabel('Time (s)')

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

-Aino



Point your RSS reader here for a feed of the latest messages in this topic.

[Privacy Policy] [Terms of Use]

© Drexel University 1994-2014. All Rights Reserved.
The Math Forum is a research and educational enterprise of the Drexel University School of Education.