Ode 45 / Ode xx: Solve time depending system of equations (2nd order)
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
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
0 Kommentare
Antworten (1)
Star Strider
am 31 Jul. 2021
Bearbeitet: Star Strider
am 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]
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,Damp,inhomo), 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,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
.
4 Kommentare
Star Strider
am 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.
.
Siehe auch
Kategorien
Mehr zu Assembly finden Sie in Help Center und File Exchange
Produkte
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!