Ode 45 / Ode xx: Solve time depending system of equations (2nd order)

Dear Matlab Forum

I´m quite new to matlab and confronted with the problem of solving a vibrational analysis (mechanics) for a car.

The linearized matrizes are the following:

Mass =  
2340           0
   0         780

Stif = 
[580000.0,                                                                                                                                                           -56000.0]
[-56000.0, 16.601100903723287828849213033975*cos(167.55160819145563938467431377491*t) - 101.84644178696417377835505444123*sin(167.55160819145563938467431377491*t) + 579200.0]

Damp =
[350.0,                                                                                                                                                                 80.0]
[ 80.0, 304.0 - 0.058599317415265314024125661602583*sin(167.55160819145563938467431377491*t) - 0.0022646900205124496311600637333048*cos(167.55160819145563938467431377491*t)]
 
inhomo =
1807.42953175056940877716695669*cos(83.775804095727819692337156887453*t) - 2326.6281369733147387054507085074*sin(83.775804095727819692337156887453*t) - 22955.400000000001455191522836685
                                - 3499.5502593343430479837509682208*cos(83.775804095727819692337156887453*t) - 3468.1908803283783962962863361566*sin(83.775804095727819692337156887453*t)

All of these matrizes are depending on the time variable t and form a system of PDEs. Therefore I can only solve the problem if t is set to 0. Otherwise i recieve an error. I tried to understand the description on the Matlab-Website where it is written, that it is possible to use a time depending Matrix M(t,y), but it wont accept my ODEs of 2nd order transformed to ODEs of 1st order.

dxdt(3:4, :) = Mass\(-inhomo - Damp*x(3:4,1) - Stif*x(1:2,1));

Does somebody knows the solution to my problem?

Here is the error description if t = 0 is commented

Unable to perform assignment because value of type 'sym' is not convertible to 'double'.

Error in testode>odefcn (line 39)
    dxdt(3:4, :) = Mass\(-inhomo - Damp*x(3:4,1) - Stif*x(1:2,1));

Error in testode>@(t,x)odefcn(t,x) (line 31)
[t,x]   = ode45(@(t,x) odefcn(t,x), tspan, [v0 a0])

Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)
  odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in testode (line 31)
[t,x]   = ode45(@(t,x) odefcn(t,x), tspan, [v0 a0])

Caused by:
    Error using symengine
    Unable to convert expression containing symbolic variables into double array. Apply 'subs' function first to substitute values for variables.
 

Here's the code I wrote.

clear, clc;

syms t
syms yC_2 yC_1 yC
syms phi_2 phi_1 phi

global Mass Stif Damp inhomo

Mass = [2340,   0;
        0,    780]
Stif = [580000,                                                                                                                                                                                                                                                                                       -56000;
        -56000, (8192*pi^2*cos((4*pi*(20*t + 4/5))/3)^2)/225 - (8192*pi^2*cos((16*pi)/15 + (80*pi*t)/3)^2)/225 - (104*pi*sin((16*pi)/5 - (160*pi*t)/3))/5 - (256*pi*sin((32*pi)/15 + (160*pi*t)/3))/15 - (1664*pi^2*cos((8*pi)/5 - (80*pi*t)/3)^2)/25 + (1664*pi^2*cos((4*pi*(20*t - 6/5))/3)^2)/25 + 579200];
Damp = [350,                                                                                                                                                                                                                                                                              80;
        80, (32*pi^2*cos((4*pi*(20*t + 4/5))/3)^2)/1125 - (32*pi^2*cos((16*pi)/15 + (80*pi*t)/3)^2)/1125 - (pi*sin((16*pi)/5 - (160*pi*t)/3))/125 - (pi*sin((32*pi)/15 + (160*pi*t)/3))/75 - (16*pi^2*cos((8*pi)/5 - (80*pi*t)/3)^2)/625 + (16*pi^2*cos((4*pi*(20*t - 6/5))/3)^2)/625 + 304];
inhomo =   [3200*sin((4*pi*(20*t + 4/5))/3) + 2600*sin((4*pi*(20*t - 6/5))/3) - 114777/5;
            3120*sin((8*pi)/5 - (80*pi*t)/3) + 2560*sin((16*pi)/15 + (80*pi*t)/3) + 2496*pi*cos((8*pi)/5 - (80*pi*t)/3) + (4096*pi*cos((16*pi)/15 + (80*pi*t)/3))/3 - (4096*pi*cos((4*pi*(20*t + 4/5))/3))/3 - 2496*pi*cos((4*pi*(20*t - 6/5))/3)];

t = 0;   % comment for error     
        
Stif = vpa(simplify(expand(subs(Stif))))
Damp = vpa(simplify(expand(subs(Damp))))
inhomo = vpa(simplify(expand(subs(inhomo))))
   
     
tspan = [0 0.1];

q0    = [0 0];
v0    = [0 0];
a0    = [0 0];

[t,x]   = ode45(@(t,x) odefcn(t,x), tspan, [v0 a0])

function dxdt = odefcn(t,x)
    global Mass Stif Damp inhomo
      
    dxdt = zeros(4,1);
    
    dxdt(1:2, :) = x(3:4,:);
    dxdt(3:4, :) = Mass\(-inhomo - Damp*x(3:4,1) - Stif*x(1:2,1));
    
end

Answer 1

[Star Strider]

The clear and clc calls are not necessary, global variables are never advisable, and the syms call is also not necessary. Create the matrices that are functions of ‘t’ as anonymous functions, and pass them as arguments to ‘odefcn’.

Try this slightly edited version of your code (plot call added) —

% clear, clc;

% syms t
% syms yC_2 yC_1 yC
% syms phi_2 phi_1 phi

% global Mass Stif Damp inhomo

Mass = [2340,   0;
        0,    780]

Mass = 2×2

2340 0 0 780

Stif_fcn = @(t) [580000,                                                                                                                                                                                                                                                                                       -56000;
        -56000, (8192*pi^2*cos((4*pi*(20*t + 4/5))/3)^2)/225 - (8192*pi^2*cos((16*pi)/15 + (80*pi*t)/3)^2)/225 - (104*pi*sin((16*pi)/5 - (160*pi*t)/3))/5 - (256*pi*sin((32*pi)/15 + (160*pi*t)/3))/15 - (1664*pi^2*cos((8*pi)/5 - (80*pi*t)/3)^2)/25 + (1664*pi^2*cos((4*pi*(20*t - 6/5))/3)^2)/25 + 579200];
Damp_fcn = @(t) [350,                                                                                                                                                                                                                                                                              80;
        80, (32*pi^2*cos((4*pi*(20*t + 4/5))/3)^2)/1125 - (32*pi^2*cos((16*pi)/15 + (80*pi*t)/3)^2)/1125 - (pi*sin((16*pi)/5 - (160*pi*t)/3))/125 - (pi*sin((32*pi)/15 + (160*pi*t)/3))/75 - (16*pi^2*cos((8*pi)/5 - (80*pi*t)/3)^2)/625 + (16*pi^2*cos((4*pi*(20*t - 6/5))/3)^2)/625 + 304];
inhomo_fcn = @(t)  [3200*sin((4*pi*(20*t + 4/5))/3) + 2600*sin((4*pi*(20*t - 6/5))/3) - 114777/5;
            3120*sin((8*pi)/5 - (80*pi*t)/3) + 2560*sin((16*pi)/15 + (80*pi*t)/3) + 2496*pi*cos((8*pi)/5 - (80*pi*t)/3) + (4096*pi*cos((16*pi)/15 + (80*pi*t)/3))/3 - (4096*pi*cos((4*pi*(20*t + 4/5))/3))/3 - 2496*pi*cos((4*pi*(20*t - 6/5))/3)];

t = 0;   % comment for error     
    
Stif = Stif_fcn(t)

Stif = 2×2

1.0e+05 * 5.8000 -0.5600 -0.5600 5.7922

Damp = Damp_fcn(t)

Damp = 2×2

350.0000 80.0000 80.0000 303.9977

inhomo = inhomo_fcn(t)

inhomo = 2×1

1.0e+04 * -2.1148 -0.3500

        
% Stif = vpa(simplify(expand(subs(Stif))))
% Damp = vpa(simplify(expand(subs(Damp))))
% inhomo = vpa(simplify(expand(subs(inhomo))))
   
     
tspan = [0 0.1];

q0    = [0 0];
v0    = [0 0];
a0    = [0 0];

[t,x]   = ode45(@(t,x) odefcn(t,x,Mass,Stif,Damp,inhomo), tspan, [v0 a0])

t = 57×1

1.0e+00 * 0 0.0000 0.0000 0.0000 0.0000 0.0001 0.0001 0.0001 0.0001 0.0003

x = 57×4

0 0 0 0 0.0000 0.0000 0.0001 0.0000 0.0000 0.0000 0.0001 0.0000 0.0000 0.0000 0.0002 0.0001 0.0000 0.0000 0.0002 0.0001 0.0000 0.0000 0.0005 0.0002 0.0000 0.0000 0.0007 0.0003 0.0000 0.0000 0.0010 0.0005 0.0000 0.0000 0.0012 0.0006 0.0000 0.0000 0.0025 0.0012

figure
plot(t, x)
grid
legend(compose('x_%d',1:size(x,2)), 'Location','best')

function dxdt = odefcn(t,x,Mass,Stif,Damp,inhomo)
%     global Mass Stif Damp inhomo
      
    dxdt = zeros(4,1);
    
    dxdt(1:2, :) = x(3:4,:);
    dxdt(3:4, :) = Mass\(-inhomo - Damp*x(3:4,1) - Stif*x(1:2,1));
    
end

.

[LeBer]

Thanks for your answer Star Strider. I implemented the tipps in the main program, but the issue regarding t still exist, since your code also only works if t = 0 is not commented.

As soon as I comment t the following happens:

% t = 0;   % comment for error 


Unrecognized function or variable 't'.

Error in testode4 (line 19)
Stif = Stif_fcn(t)

I´m not sure if I understood the ode functionality completly. Is it possible to solve ode's with respect to a certain variable, e.g. t as in this example?

[Star Strider]

I am not certain what your equation is.

One option is to use the pdepe function to integrate it. I have used it a few times, although not recently, so I would have to know a bit more about the function you want to integrate in order to possibly help you to use pdepe to solve it. (I am even less familiar with the Partial Differential Equation Toolbox since I have not used it for at least the last decade. I will do my best to help you with it if necessary, however it would take me some time to understand in some detail in order to be minimally proficient in how to use it.)

Otherwise, there are many ways to define the time variable and use it with your system using the ODE integration functions. It depends what you want to do.

.

[LeBer]

Hmm, I´ll try to explain my problem. Thanks for the pdepe offer, but I think the ode-function could solve it.

1) My two ode functions are coupled, linear (after I linearized them) and time depending

2) the matrizes (Mass, Damp, Stif, inhomo) are derived from the linearized problem

3) in the matrizes there is still a time dependency. So the 2 odes change with time.

I marked the first time derivative by _1 (= dot) and the 2nd time derivative by_2 ( = dotdot)

Mass * [x1_2;x2_2] + Damp(t) * [x1_1;x2_1] + Stif(t) * [x1; x2] = inhomo(t)

Now I need to solve this equation w.r.t. x1, x2 and t, but it only works if t is set to a specific number (e.g 0 as seen above). The odeXX function should vary t while solving the ode for the specific time steps, like the docs say ( "All MATLAB® ODE solvers can solve systems of equations of the form y′=f(t,y), or problems that involve a mass matrix, M(t,y)y′=f(t,y)")

% vary t and solve with odeXX
Mass * [x1_2;x2_2] + Damp(0.01) * [x1_1;x2_1] + Stif(0.01) * [x1; x2] = inhomo(0.01)
...
Mass * [x1_2;x2_2] + Damp(1.5) * [x1_1;x2_1] + Stif(1.5) * [x1; x2] = inhomo(1.5)

So my question is: How do I make odeXX change the time-variable t while solving?

[Star Strider]

If I understand correctly what you want to do, this is likely straightforward, described in ODE with Time-Dependent Terms so well-established. However, that may not be necessary here.

Since they are all functions of time, evaluate them as such within ‘odefcn’.

Try this —

% global Mass Stif Damp inhomo

Mass = [2340,   0;
        0,    780];
Stif_fcn = @(t) [580000,                                                                                                                                                                                                                                                                                       -56000;
        -56000, (8192*pi^2*cos((4*pi*(20*t + 4/5))/3)^2)/225 - (8192*pi^2*cos((16*pi)/15 + (80*pi*t)/3)^2)/225 - (104*pi*sin((16*pi)/5 - (160*pi*t)/3))/5 - (256*pi*sin((32*pi)/15 + (160*pi*t)/3))/15 - (1664*pi^2*cos((8*pi)/5 - (80*pi*t)/3)^2)/25 + (1664*pi^2*cos((4*pi*(20*t - 6/5))/3)^2)/25 + 579200];
Damp_fcn = @(t) [350,                                                                                                                                                                                                                                                                              80;
        80, (32*pi^2*cos((4*pi*(20*t + 4/5))/3)^2)/1125 - (32*pi^2*cos((16*pi)/15 + (80*pi*t)/3)^2)/1125 - (pi*sin((16*pi)/5 - (160*pi*t)/3))/125 - (pi*sin((32*pi)/15 + (160*pi*t)/3))/75 - (16*pi^2*cos((8*pi)/5 - (80*pi*t)/3)^2)/625 + (16*pi^2*cos((4*pi*(20*t - 6/5))/3)^2)/625 + 304];
inhomo_fcn = @(t)  [3200*sin((4*pi*(20*t + 4/5))/3) + 2600*sin((4*pi*(20*t - 6/5))/3) - 114777/5;
            3120*sin((8*pi)/5 - (80*pi*t)/3) + 2560*sin((16*pi)/15 + (80*pi*t)/3) + 2496*pi*cos((8*pi)/5 - (80*pi*t)/3) + (4096*pi*cos((16*pi)/15 + (80*pi*t)/3))/3 - (4096*pi*cos((4*pi*(20*t + 4/5))/3))/3 - 2496*pi*cos((4*pi*(20*t - 6/5))/3)];

% t = 0;   % comment for error     
%     
% Stif = Stif_fcn(t)
% Damp = Damp_fcn(t)
% inhomo = inhomo_fcn(t)
        
% Stif = vpa(simplify(expand(subs(Stif))))
% Damp = vpa(simplify(expand(subs(Damp))))
% inhomo = vpa(simplify(expand(subs(inhomo))))
   
     
tspan = [0 0.1];

q0    = [0 0];
v0    = [0 0];
a0    = [0 0];

[t,x]   = ode45(@(t,x) odefcn(t,x,Mass,Stif_fcn,Damp_fcn,inhomo_fcn), tspan, [v0 a0]);


figure
plot(t, x)
grid
legend(compose('x_%d',1:size(x,2)), 'Location','best')

function dxdt = odefcn(t,x,Mass,Stif_fcn,Damp_fcn,inhomo_fcn)
%     global Mass Stif Damp inhomo
      
    dxdt = zeros(4,1);
    
    dxdt(1:2, :) = x(3:4,:);
    dxdt(3:4, :) = Mass\(-inhomo_fcn(t) - Damp_fcn(t)*x(3:4,1) - Stif_fcn(t)*x(1:2,1));
    
end

Note — The functions themselves are now being passed as arguments, and evaluated at each time step.

.