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

6 views (last 30 days)
LeBer on 31 Jul 2021
Commented: Star Strider on 3 Aug 2021
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

Star Strider on 31 Jul 2021
Edited: Star Strider on 31 Jul 2021
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
.
Star Strider on 3 Aug 2021
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.
.