Solving Differential Equations With

Hi, I'm new in Matlab and I have a problem with ODE45 comand.

Basically, I'm trying to solve the second order differential equation of a spring (ke) -mass (me) -damper (ce) system. With [0; 0] initial conditions. It follows the code:

[t,x] = ode45(@SMD,[0 Tmax],[0; 0]);
   plot(t,x(:,1))

SMD is defined in a function as:

function dxdt = SMD(t,x,me,ce,ke,Ft)
fe=2;
Qsie=0.05;
me=32000;
ce=2*Qsie*(2*pi*fe)*me;
ke=me*(2*pi*fe)^2;
dxdt = [x(2); (-ce*x(2)-ke*x(1)+Ft)/me];

But Ft is a force vector (varying on time) calculated previously. How can I import this vector to dxdt function?

Is there any restriction in relation to the Length of vector Ft? I asked this question because I dont know the time increment of the integration when I use ODE45. The vectors of the analysis may be different Lengths (Ft and x).

Other question: there is a form of the variables me, ce and ke be imported from the main program? When I try compile the code without indicate the values of them in the function dxdt, an error was reported.

Thanks very much!

[Torsten]

What is the time vector corresponding to Ft ?

Answer 1

[Sam Chak]

Hi @Igor Braz Gonzaga

I add a new answer because the pieces of new info were scattered here and there. You should have edited the original post and included the new info about the stochastic load. Anyhow, I've added the interpolation (as initially suggested by @Torsten and @Walter Roberson) in the odefcn to get the stochastically continuous load.

function dydt = odefcn(t, y, tF, Fft)
    dydt = zeros(2,1);
        fe = 2;
        Omega = 2*pi*fe;                % undamped natural frequency
        Qsie = 0.05;                    % damping ratio
        me = 32000;                     % “generalized mass” coefficient
        ce = 2*Qsie*(Omega)*me;         % damping coefficient 
        ke = me*(Omega)^2;              % stiffness coefficient 
        A = [0 1; -ke/me -ce/me];       % state matrix from ODE
        B = [0; 1/me];                  % input matrix from ODE
        p = [-1+1i -1-1i];              % eigenvalues from desired x" + 2*x' + 2*x = 0
        K = place(A, B, p);             % control gains
        Ft1 = - K(1)*y(1) - K(2)*y(2);  % Case #1: control force
        Ft2 = 0;                        % Case #2: no force (free response)
        Ft3 = me*cos(Omega*t);          % Case #3: disturbing sinusoidal force
        Ft4 = interp1(tF, Fft, t);      % Interpolate the stochastic data set
    dydt(1) = y(2);
    dydt(2) = (- ce*y(2) - ke*y(1) + Ft4)/me;
end

The stochastic load is generated from the teste_exluir m-file, but some variables are undefined. So, I have assigned some values to the number of pedestrians, Nped and L (because you never specify what it is), and renamed the time and force variables to tF and Fft, respectively. Rest assured that the stochastic load formula remains unchanged. I have added a few lines at the end of the script and ran for 3 tests. Please review if the results are expected or acceptable for your study.

Fft = Ft(:,1);      % Force induced by the walking pedestrians
figure(1)
plot(tF, Fft)

tspan = [0 10];
y0 = [0 0];
[t, y] = ode45(@(t, y) odefcn(t, y, tF, Fft), tspan, y0);
figure(2)
plot(t, y)

Results:

[Walter Roberson]

You cannot use the default interp1() method for this purpose: you need to use something like 'spline' or 'cubic' (not sure about 'pchip')

[Torsten]

Further, assuming the integration scheme could handle such discontinuous inputs (which it can't), this is one possible realization of your stochastic process for Ft. The next realization will yield a different response. So you would have to run your equations thousands of times to get the distribution (not the values) of your unknowns y(1) and y(2).

That's where stochastic differential equations and the associated methods come into play.

[Sam Chak]

Yes, @Igor Braz Gonzaga literally has to manually run the stochastic differential equations (SDE) many times with this approach, if he wants to perform the Monte Carlo Simulation. Unfortunately, I don't have the Financial Toolbox to perform SDE. I have checked the pricing here, if he wishes to get the license.

https://www.mathworks.com/pricing-licensing.html?prodcode=FI

By the way, is there a workaround for this problem that @Walter Roberson and @Torsten can suggest to him so that he can automate the simulations and storing the data for a range of values for Nped and L in the teste_exluir script using the Standard MATLAB? Thank you for considering this request.

Answer 2

[Torsten]

fe=2;
Qsie=0.05;
me=32000;
ce=2*Qsie*(2*pi*fe)*me;
ke=me*(2*pi*fe)^2;
[t,x] = ode45(@(t,x)SMD(t,x,me,ce,ke,tt,Ft),[0 Tmax],[0; 0]);

function dxdt = SMD(t,x,me,ce,ke,tt,ft)
  Ft = interp1(tt,ft,t);
  dxdt = [x(2); (-ce*x(2)-ke*x(1)+Ft)/me];
end

where tt is the time vector corresponding to Ft.

[Igor Braz Gonzaga]

Hi Torsten, thanks a lot for your comment. But an error was reported using your code:

Mensage error: " Error in ode45 (line 115) odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin) "

My code follows:

Qsie=0.05;
me=32000;
ce=2*Qsie*(2*pi*fe)*me;
ke=me*(2*pi*fe)^2;

[temp,x] = ode45(@(temp,x)SMD(temp,x,me,ce,ke,t,Ft),[0 Tmax],[0; 0]);

where t is a vector t=0:1/(100*fe):Tmax; fe=2 (integer) and Ft is a vector of dimension length(t).

The SMD function:

function dxdt = SMD(temp,x,me,ce,ke,t,ft)
Ft = interp1(t,ft,temp);
dxdt = [x(2); (-ce*x(2)-ke*x(1)+Ft)/me];
end

Please, If you could help me. I will be thankful.

[Torsten]

Except that "fe" and "Tmax" do not appear in your list of parameter definitions, I don't see an error in your code.

[Igor Braz Gonzaga]

Hi Torsten, I think that the error is on variable Tmax. When I define this variable with a value of 20 seconds, the code run perfectly. But if I change this value, for example, 30 seconds, an mensage error is reported by Matlab.

%Structure date
fe=2;
Qsie=0.05;
me=32000;
ce=2*Qsie*(2*pi*fe)*me;
ke=me*(2*pi*fe)^2;
%Time date
Tmax=20; %change the value of this variable reports a problem
t=0:1/50:Tmax;
Ft=50*750*sin(tt*5);

[temp,x] = ode45(@(temp,x)solve_(temp,x,me,ce,ke,t,Ft),[0 50],[0; 0]);
plot(temp,x(:,2))

function dxdt = solve_(temp,x,me,ce,ke,t,ft)
Ft = interp1(t,ft,temp);
dxdt = [x(2); (-ce*x(2)-ke*x(1)+Ft)/me];
end

When Tmax is 20 s, Matlab plots the following graph (I did not understand this result, the variable Temp was set from 0 to 50, does not this chart have to plot the time up to 50 s? ):

The following mensage error appears when I change Tmax (change from 20s):

If you know what is the problem, please, help me.

Thanks.

[Walter Roberson]

You do not define tt for use in Ft, but if you replace tt with t then I see no error message.

However if you switch to ode15s you would get an error. That is because you are doing linear interpolation for Ft, but linear interpolation breaks the mathematics of the Runge-Kutta methods, which require that the second derivatives of the functions are continuous.

%Structure date
fe=2;
Qsie=0.05;
me=32000;
ce=2*Qsie*(2*pi*fe)*me;
ke=me*(2*pi*fe)^2;
%Time date
Tmax=30; %change the value of this variable reports a problem
t=0:1/50:Tmax;
Ft=50*750*sin(t*5);

[temp,x] = ode45(@(temp,x)solve_(temp,x,me,ce,ke,t,Ft),[0 50],[0; 0]);
plot(temp,x(:,2))

function dxdt = solve_(temp,x,me,ce,ke,t,ft)
Ft = interp1(t,ft,temp);
dxdt = [x(2); (-ce*x(2)-ke*x(1)+Ft)/me];
end

[Igor Braz Gonzaga]

Thanks Walter Roberson, Its done!

But now I have some other problem, I change Ft to be a stochastic load. The dimension of Ft continues to be the same (length(t)). The following graph shows Ft in time.

I dont know what is happen, but with this change of Ft, an error is reported:

If you know how to fix this problem, please, help me.

Thanks a lot.

[Torsten]

This is a general remark, not directly related to the error you receive:

The function Ft cannot be stochastic if you use an ODE solver from MATLAB's solver suite.

The inputs for the solvers are assumed at least continuous. So the result of using the Ft curve from above will either be that the solver stops integration or that the results will be arbitrary.

For the error message, please include your complete code.

[Igor Braz Gonzaga]

Hi Torsten, In this case, for a random force (discrete in time), is there any command in Matlab to solve this ODE?

[Walter Roberson]

@Igor Braz Gonzaga

None of the ode*() functions support discontinuities in the equations such as might be due to random forces or impulses.

The Finance toolbox has some support for Stochastic Differential Equations; https://www.mathworks.com/help/finance/stochastic-differential-equation-sde-models.html

For the ode*() routines, you would need to stop integrating at every point that randomness occurs, and resume integrating again.

Answer 3

[Sam Chak]

Hi @Igor Braz Gonzaga

If this is an exercise in the Mechanics class, usually will be given. Anyhow, there are generally 3 cases for the force vector, :

• Case 1: feedback control force;
• Case 2: no force input (free response);
• Case 3: external disturbing force (can be a sinusoidal signal or a constant force).

You can investigate the system response by selecting the type of force (Ft1, Ft2, or Ft3) on the “dydt(2)” line.

function dydt = odefcn(t, y)
    dydt = zeros(2,1);
        fe = 2;
        Omega = 2*pi*fe;                % undamped natural frequency
        Qsie = 0.05;                    % damping ratio
        me = 32000;                     % mass coefficient
        ce = 2*Qsie*(Omega)*me;         % damping coefficient 
        ke = me*(Omega)^2;              % stiffness coefficient 
        A = [0 1; -ke/me -ce/me];       % state matrix derived from ODE (for Case #1 only)
        B = [0; 1/me];                  % input matrix derived from ODE (for Case #1 only)
        p = [-1+1i -1-1i];              % eigenvalues from desired x" + 2*x' + 2*x = 0 (for Case #1 only)
        K = place(A, B, p);             % control gains (for Case #1 only)
        Ft1 = - K(1)*y(1) - K(2)*y(2);  % Case #1: control force
        Ft2 = 0;                        % Case #2: no force (free response)
        Ft3 = me*cos(Omega*t);          % Case #3: disturbing sinusoidal force
    dydt(1) = y(2);
    dydt(2) = (- ce*y(2) - ke*y(1) + Ft1)/me;
end

You can also set the simulation time, tspan, and the initial condition, y0, depending on your original problem. (, ) is the equilibrium point for Cases #1 & #2. If this is a control project, you will need to determine the eigenvalues (a.k.a poles) from the desired system response and then the control gains can be designed using the pole placement technique.

https://www.mathworks.com/help/control/ref/place.html

tspan = linspace(0, 20, 2001)';
y0 = [0.1; 0];
[t, y] = ode45(@(t, y) odefcn(t, y), tspan, y0);
plot(t, y)

Results:

[Igor Braz Gonzaga]

Hi Sam, thanks a lot for your comment. My case is The 3rd. But in my case, the force is not a harmonic load, it is a random load (discrete). Do you know how to solve this Ode?

[Sam Chak]

Hi @Igor Braz Gonzaga

I didn't notice that you posted in the comment section. But this looks like a harmonic load. Can you reconfirm the math equation for ? If it is a critical information, you can actually edit your Question and put the new info there.

Else, please clarify if you are free to assign as any discrete-time random load, or if can be converted to a Stairstep signal.

https://www.mathworks.com/help/matlab/ref/stairs.html

[Igor Braz Gonzaga]

Hi @Sam Chak,

Initially I just tested a simple harmonic load Ft=(50)*750*sin(5t), but I have to solve a problem in which the external disturbing force is not a deterministic load (such as the sinusoidal harmonic load that I used).

I've already implemented my random disturbing force that is given by the following code:

for i=1:Nped
    for j=1:lenght(t)
        for k=1:NH
            Ft=G(i)+G(i)*alfa(k)*sin(2*pi*i*fp(i)*t(j)+phase(k))
        end
    end 
end

where: G(i) is a vector of the weight of the N-pedestrians (Nped); fp(i) is a vector of the walking frequency of the N-pedestrians; alfa(k) and phase(k) are vectors of constant values. The weight G(i) and walking frequency fp(i) of each pedestrian was randomly generated for each person.

The following graph shows an example of Ft (there are other steps of calculus not showed before, I dont put it on the code to simplify):

So, when I try to run my code with an harmonic load it works well, but when I input this random Ft on the ODE45 command, it does not work =/

Could you help me?

Thanks so much!!

[Torsten]

Are you sure about your loop from above ? It will only give you one single value for Ft since Ft is a scalar.

[Igor Braz Gonzaga]

I'm sorry Torsten. There is an error. Here I paste the correct code:

for i=1:Nped
    for j=1:length(t)
        for k=1:NH       
            Findiv(i,j)=(G(i)+alfadin(k)*G(i)*sin(2*pi*i*fp(i)*(t(j)-Tentr(i))-phase(k)));
        end
        Ft(j)=Ft(j)+Findiv(i,j);
    end
end

[Sam Chak]

Hi @Igor Braz Gonzaga

Thanks for your reply. I guess you are trying to model the discrete oscillation that is generated by random number of pedestrians walking over a rigid body structure. Since you can freely model it, can you show the desired equation in LaTeX format?

[Igor Braz Gonzaga]

Of Course @Sam Chak.

where is the modal shape of the structure in relation to the instantaneous position of the person on the structure.

The following graph shows an example of :

[Torsten]

Best you include the data directly as a .mat -file or a .txt file.

But as I said:

You can't use stochastic inputs in the ODE solvers.

Fit a bell-shaped curve to your data - this might work then.

[Igor Braz Gonzaga]

@Torsten, follow the files attachment.

So, if I cant use ODEsolvers for stochastic inputs, there is another command on Matlab to do that?

[Sam Chak]

Hi @Igor Braz Gonzaga

Thanks for the input files. As suggested by @Torsten, I have made an approximation of the Ft using the bell curve or normal distribution function with added noise. Is this good enough for the vibration analysis of the structure? It is also possible to make a skewed distribution shape, if needed. I think using the Gamma distribution function.

t = [0:0.02:15]';
Amp = 3.5e4;
mean = 7.5
stddev = 1.5
x = Amp*stddev*sqrt(2*pi)*normpdf(t, mean, stddev);
y = awgn(x, 30, 'measured');
plot(t, [x y])

This way, you should be able to run ODE45.

[Igor Braz Gonzaga]

Hi Sam Chak, I didn't understand. The input load that I am generating is a moving force of random pedestrians walking on a rigid structure. For that, the parameters of weight (G) and walking frequency (fp) of the pedestrians should be random assigned.

Nevertheless, just this parameters follows normal distribution (they are obtained from normal distribution). The force Ft generated by these random pedestrians do not follows a normal distribution shape.

The following graph shows an example of Ft for Nped=20.

[Torsten]

As pointed out several times, you cannot feed the differential equation with stochastic input.

So either you accept a way to smoothen the Ft curves or we should quit discussion here.

[Walter Roberson]

As I posted earlier,

The Finance toolbox has some support for Stochastic Differential Equations; https://www.mathworks.com/help/finance/stochastic-differential-equation-sde-models.html

If you do not use that, then you need to stop integration each time you have a new random event.

Answer 4

[Igor Braz Gonzaga]

Hi @Sam Chak, thanks a lot!! It works!!

Yes, My ideia is do a loop with a series of simulations (Monte Carlo) to attain the responses of the structure under stochastic load.

But I have another question, more simple, I guess.

If I want the acceleration of the structure, what should I do?

I try to derive y(:,2) vector, using the following code:

acel=diff(y(:,2));
t_acel=(tt(2:end)+tt(1:(end-1)))/2;
plot(t_acel,acel)

but the graph do not correspond with the expected values , that is, analytical solution (I try to use a more simple example, with a deterministic input load).

The graphs of y(:,2) (velocity of the structure) and y(:,1) (displacement of the structure) works well, but the acceleration does not.

Thanks

[Walter Roberson]

acel=diff(y(:,2));

That assumes that the items are all the same time apart.

You should use gradient(y(:,2), tt)