ode45 forcing function

Question

Ganesh 2011-8-18

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/13959-ode45-forcing-function

编辑： mbvoyager 2018-7-26

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 个评论
显示 -1更早的评论隐藏 -1更早的评论

Jan 2011-8-18

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 2018-7-26

3
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/13959-ode45-forcing-function#answer_330485

编辑：mbvoyager 2018-7-26

在 MATLAB Online 中打开

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