ode45 forcing function

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

[Jan]

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.

Answer 1

[mbvoyager]

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.

Answer 2

[Floris]

I think this is what you need:

http://www.mathworks.se/help/techdoc/ref/ode23.html Example 3

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...