MATLAB Answers

ode45 forcing function

119 views (last 30 days)
Ganesh
Ganesh on 18 Aug 2011
Edited: mbvoyager on 26 Jul 2018
I am solving a sdof system with the forcing function as a random process.I have the forcing function as a vector.I am having difficulties to solve the equation with this.Can anyone please provide me some help on how i can input the forcing vector and thus solve the equation using ode45?
Thanks in advance.
Ganesh

  1 Comment

Jan
Jan on 18 Aug 2011
What is a "forcing function"? And what do you expect, when you integrate a system based on a random process. ODE45 needs a smooth function, otherwise the result is random also.

Sign in to comment.

Answers (2)

mbvoyager
mbvoyager on 26 Jul 2018
Edited: mbvoyager on 26 Jul 2018
Recently I had to accomplish the same task.
The "issue" is that ode45 of course changes the time step size due to the stability criteria. It implements a adaptive time step size.
This means that you need values from your stochastic process according to this adaptive time steps.
An easy way is to interpolate between the discretized values of your stochastic process and the demanded time of the ode45 function.
For demonstration I used an example from the Matlab ode45 documention, and applied a brownian motion to this model that can be seen as an artificial parameter.
intB
So to combine those two models, this is what I've done:
%%Model Parameters
%Start and End Time
tspan = [0 5];
%Time discretization for Stochastic Process
dt=0.001;
%Time Domain for Stochastic Process
tt=0:dt:tspan(end);
%Boundary Condition
y0 = 0;
%%Stochastic Process
%(http://www-math.bgsu.edu/~zirbel/sde/matlab/brownian.m)
%Drift coefficient
b=1;
%Diffusion coefficient
sigma=10;
%Model of a Brownian Motion
W = [0; cumsum(randn(length(tt)-1,1))]/sqrt(length(tt)-1);
W = W*sqrt(tt(end));
B = b*tt' + sigma*W;
%Plot the Stochastic Process
figure
plot(tt,B)
hold on
plot(tt,b*tt,':');
title([int2str(length(tt)) '-step version of Brownian motion and its mean'])
xlabel(['Drift ' num2str(b) ', diffusion coefficient ' num2str(sigma)])
%%Solving 2nd order van der Pol
%(https://de.mathworks.com/help/matlab/ref/ode45.html?searchHighlight=ode45&s_tid=doc_srchtitle)
%Solving procedure of the function using FDM of ode45 including a
% stochastic process
[t,y] = ode45(@(t,y) vdp1(t,y,tt,B),tspan,[2; 0]);
figure
plot(t,y(:,1),'-o',t,y(:,2),'-o')
title('Solution of van der Pol Equation (\mu = 1) with ODE45');
xlabel('Time t');
ylabel('Solution y');
legend('y_1','y_2')
Now this is the crucial part, for every call of ode45 to the function vdp1, a value from the stochastic process needs to be interpolated.
function dydt = vdp1(t,y, tt, B)
%VDP1 Evaluate the van der Pol ODEs for mu = 1
%
% See also ODE113, ODE23, ODE45.
% Jacek Kierzenka and Lawrence F. Shampine
%Interpolation of Stochastic Process
intB = interp1(tt,B,t);
dydt = [y(2); (1-y(1)^2)*y(2)-y(1)-intB];
end
This is the solution of the ordinary 2nd order van der Pol equation (seen in the matlab documentation example):
And this is the solution where a stochastic process has induced some changes in every time step:
Copyright 1984-2014 The MathWorks, Inc.

  0 Comments

Sign in to comment.


Floris
Floris on 6 Sep 2011
I think this is what you need:
where the ODE receives a time-dependent function, defined by the user. I think this is your "forcing function".
But a "random process" will give random results as well. Plus, I'd be concerned about the stability of the solution...

  0 Comments

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by