Problem in solving matrix diiferential equation using ode45?
Show older comments
Where is the problem in my code I could not solve a matric differential equations ?
M=[7.92982501567794e-08 7.16408455108792e-12 -1.30449341090077e-11 -1.67007458454893e-12 -1.69491338826731e-12 2.04671335711028e-12 -5.39085143946466e-12 -6.86316937563370e-12 1.13976593936898e-12 8.73221272415655e-12 7.29148797499054e-12
7.16408455108792e-12 7.92999948060020e-08 -9.01780276368023e-09 -4.19265297003500e-11 -1.21828612766692e-11 -2.89846125284266e-12 -2.48017467491464e-12 2.65328942481270e-12 -5.79570472068573e-13 6.17511954582720e-12 -4.72437536818714e-13
-1.30449341090077e-11 -9.01780276368023e-09 7.92999125348644e-08 4.00936687032002e-11 2.04283104742536e-11 1.24820735380079e-12 9.30760973993229e-12 -7.50079609806551e-12 1.22865516464498e-12 8.06254006030017e-12 5.56252493353000e-12
-1.67007458454893e-12 -4.19265297003500e-11 4.00936687032002e-11 7.92994152491339e-08 -2.20189255359344e-12 -1.19756353847227e-11 9.16324192448759e-12 -8.26074802470566e-12 -2.45588226225298e-12 4.46842214440270e-13 3.39106218157172e-12
-1.69491338826731e-12 -1.21828612766692e-11 2.04283104742536e-11 -2.20189255359344e-12 7.92999951227102e-08 1.69566385897829e-08 -1.13259218230788e-11 -1.06306934205658e-11 5.69865748002221e-13 -1.33330158047386e-11 -2.77602627891578e-12
2.04671335711028e-12 -2.89846125284266e-12 1.24820735380079e-12 -1.19756353847227e-11 1.69566385897829e-08 7.92999831052483e-08 1.16054950164138e-12 5.22206051556241e-12 -7.72954490091433e-12 -2.47645215082267e-12 -3.12393940250136e-12
-5.39085143946466e-12 -2.48017467491464e-12 9.30760973993229e-12 9.16324192448759e-12 -1.13259218230788e-11 1.16054950164138e-12 7.92999555250990e-08 -6.52483339571384e-10 -1.19278810548402e-12 5.25758800878068e-12 6.39317043900566e-12
-6.86316937563370e-12 2.65328942481270e-12 -7.50079609806551e-12 -8.26074802470566e-12 -1.06306934205658e-11 5.22206051556241e-12 -6.52483339571384e-10 7.92994823357396e-08 3.66837619211617e-12 -2.15281354934899e-11 -2.06399667934130e-13
1.13976593936898e-12 -5.79570472068573e-13 1.22865516464498e-12 -2.45588226225298e-12 5.69865748002221e-13 -7.72954490091433e-12 -1.19278810548402e-12 3.66837619211617e-12 7.92997780363045e-08 -3.20478613318390e-09 -9.20728674348331e-13
8.73221272415655e-12 6.17511954582720e-12 8.06254006030017e-12 4.46842214440270e-13 -1.33330158047386e-11 -2.47645215082267e-12 5.25758800878068e-12 -2.15281354934899e-11 -3.20478613318390e-09 7.92999432448874e-08 4.27648531379977e-12
7.29148797499054e-12 -4.72437536818714e-13 5.56252493353000e-12 3.39106218157172e-12 -2.77602627891578e-12 -3.12393940250136e-12 6.39317043900566e-12 -2.06399667934130e-13 -9.20728674348331e-13 4.27648531379977e-12 7.92999522395012e-08];
C=[-3.42898926854769e-06 -3.09853079838774e-10 5.63680179774281e-10 7.18316123505485e-11 7.32043974513657e-11 -8.83641052359976e-11 2.32608057501850e-10 2.96893262006781e-10 -4.94834376380811e-11 -3.77414155634669e-10 -3.15280315155373e-10
-3.09853079838774e-10 -3.42822247644157e-06 3.89849126448742e-07 1.81237812695370e-09 5.26664903076789e-10 1.25368640746015e-10 1.07176084927423e-10 -1.15093176288670e-10 2.53119614277008e-11 -2.66878262734472e-10 2.05557932408960e-11
5.63680179774281e-10 3.89849126448742e-07 -3.42821883018528e-06 -1.73315988046962e-09 -8.82865166143735e-10 -5.37877827811127e-11 -4.02403300820132e-10 3.30719950046673e-10 -5.67461808821010e-11 -3.48413594745865e-10 -2.40423082615653e-10
7.18316123505485e-11 1.81237812695370e-09 -1.73315988046962e-09 -3.42735817917801e-06 9.52596198600402e-11 5.17542525850661e-10 -3.95943650551844e-10 3.56956034992411e-10 1.06211785588847e-10 -1.94110737204741e-11 -1.46516691611563e-10
7.32043974513657e-11 5.26664903076789e-10 -8.82865166143735e-10 9.52596198600402e-11 -3.42682637747221e-06 -7.32754904764228e-07 4.89474123375667e-10 4.59411557481755e-10 -2.46702162282727e-11 5.75962888883721e-10 1.19725151113121e-10
-8.83641052359976e-11 1.25368640746015e-10 -5.37877827811127e-11 5.17542525850661e-10 -7.32754904764228e-07 -3.42682555369805e-06 -5.00861421182035e-11 -2.25669451568986e-10 3.33860804725077e-10 1.07182548547866e-10 1.34689134354363e-10
2.32608057501850e-10 1.07176084927423e-10 -4.02403300820132e-10 -3.95943650551844e-10 4.89474123375667e-10 -5.00861421182035e-11 -3.42599149801318e-06 2.81892081456211e-08 5.16487620733699e-11 -2.26966026963761e-10 -2.76288740535447e-10
2.96893262006781e-10 -1.15093176288670e-10 3.30719950046673e-10 3.56956034992411e-10 4.59411557481755e-10 -2.25669451568986e-10 2.81892081456211e-08 -3.42597017380682e-06 -1.50581477924180e-10 9.29819231192795e-10 9.05872401930873e-12
-4.94834376380811e-11 2.53119614277008e-11 -5.67461808821010e-11 1.06211785588847e-10 -2.46702162282727e-11 3.33860804725077e-10 5.16487620733699e-11 -1.50581477924180e-10 -3.42487796344742e-06 1.38411425482107e-07 3.97409661925509e-11
-3.77414155634669e-10 -2.66878262734472e-10 -3.48413594745865e-10 -1.94110737204741e-11 5.75962888883721e-10 1.07182548547866e-10 -2.26966026963761e-10 9.29819231192795e-10 1.38411425482107e-07 -3.42488187253882e-06 -1.84679726986698e-10
-3.15280315155373e-10 2.05557932408960e-11 -2.40423082615653e-10 -1.46516691611563e-10 1.19725151113121e-10 1.34689134354363e-10 -2.76288740535447e-10 9.05872401930873e-12 3.97409661925509e-11 -1.84679726986698e-10 -3.42460989035060e-06];
K=[0.763518187446834 -0.123471983781582 2.12644497110045 -0.00159184932062271 0.00286665123297165 0.00326795662599221 -0.00565823204629348 -0.0111247074471964 -0.852444758506581 0.0170129698870622 0.000567993613223703
0.122986250969018 4.76907928529478 -0.543239600079982 2.12608568267782 0.220268589553943 0.0633458824464058 -0.00720500630476635 5.67455459314993e-05 0.000218541703671969 0.00504540957654056 0.00297905359335030
-2.13057558844554 -0.544693846047297 4.79735395984435 -0.116535843976208 -3.80931536422324 -1.08160690561449 0.00135336667261543 0.00474788044688154 0.00305299632669883 0.00159367036253352 -0.0179988437727940
-0.00372848570465477 -2.13382763678292 0.123763968250996 12.3142831297151 -0.000137283537755157 -0.00228005296293815 3.81673953925391 -0.0591314258833513 0.00127018393509301 -0.00796895715491763 0.00234941752038280
-0.00163737934937908 -0.222001483568310 3.81734112223475 0.00173003724327611 19.0472196037760 4.07106712942156 -0.00255853650434496 -0.125116265764028 5.46695676693670 -0.0993709375445011 0.00703634900851892
0.00258298057052472 -0.0623605607181842 1.08529817957274 -0.00267072015732558 4.07294423066431 19.0574340781544 0.0147753791569810 2.05513925508344 1.52210287227782 -0.0237523853866786 0.00922809527607430
-0.00468535654428222 -0.00626853319266933 -0.00115469998983706 -3.81119446621267 -0.00155959499500691 -0.0121660190928631 32.1310954783617 -0.262056094807989 0.000541439656752278 0.00464738145837403 -0.0255209825354020
0.0134976568920513 -0.000376978112019705 0.00263963824092268 0.0542232932656596 0.120732014828162 -2.05279502259563 -0.263280871097288 32.3912983853784 0.00480792646303815 -0.0114157922244635 -3.80017279650533
0.847260618562458 0.000667255021819536 0.00160008219803174 0.000312861305428450 -5.46833245082724 -1.53174909321691 0.00133283528195683 0.00562883566011632 55.0224400988869 -2.22476909528779 0.00148357487599929
-0.0121173666954115 0.00477027008481583 0.00127363569963892 0.00477817239843339 0.0854464228914320 0.0287153890552796 0.00377260121504029 -0.0113582049525324 -2.22603066376352 55.0972049682465 0.00355018925501934
0.000736959774981387 0.00110272425151785 0.0202945683248242 0.000494668031994278 0.00591913812055584 0.00913408539493314 0.0270935852755023 3.80595460544076 -0.00385777038321837 0.00320054596414318 61.4845746182663];
F=[-0.00099996606222035952475160617808552
-0.0000002460836756585278896861534876451
-0.00000060555830113300843168099357979169
-0.00000023236067999397214860565123562071
-0.00092724784232042312137180967502376
-0.0012336739219426489876333403630418
0.00000015739906629771823481279895519369
0.000000049705400083053920707370130069998
0.0000000074445232833242759638623669833003
0.0000001051612538093116976310734137023
-0.00099218997013689166348230850119706];
syms tau_1(t) tau_2(t) tau_3(t) tau_4(t) tau_5(t) tau_6(t) tau_7(t) tau_8(t) tau_9(t) tau_10(t) tau_11(t)
v = transpose([tau_1(t) tau_2(t) tau_3(t) tau_4(t) tau_5(t) tau_6(t) tau_7(t) tau_8(t) tau_9(t) tau_10(t) tau_11(t)]);
odes = M * diff(diff(v,t)) + C * diff(v,t) + K * v(t) - F ;
Y = matlabFunction(odes,'Vars',v(t));
Y0 = [0;0];
tspan = [0, 1];
[t, y] = ode45(@Y, tspan, Y0);
plot(t, y);
1 Comment
Dyuman Joshi
on 4 Jan 2024
Edited: Dyuman Joshi
on 4 Jan 2024
Your code will not work, even if it is (incorrectly) corrected in the given answer, because the symbolic functions tau_1 to tau_11 are not explicitly defined.
Accepted Answer
Hi @Mehdi
Since you already have the nicely-defined matrices, you might as well put the ODEs in the State-Space model and solve the system using the ode45 solver. However, the responses indicate that the system is unstable. It's advisable to double-check the values in the matrices.
tspan = [0, 3e-3];
v0 = zeros(22, 1); % Use your actual initial values
[t, v] = ode45(@odefun, tspan, v0);
plot(t, v(:, 1:11)), grid on % Plot the first 11 states
xlabel t, ylabel \bf{v}_{1:11}
%% State-space Model
function dvdt = odefun(t, v)
%% Mass matrix
M=[7.92982501567794e-08 7.16408455108792e-12 -1.30449341090077e-11 -1.67007458454893e-12 -1.69491338826731e-12 2.04671335711028e-12 -5.39085143946466e-12 -6.86316937563370e-12 1.13976593936898e-12 8.73221272415655e-12 7.29148797499054e-12
7.16408455108792e-12 7.92999948060020e-08 -9.01780276368023e-09 -4.19265297003500e-11 -1.21828612766692e-11 -2.89846125284266e-12 -2.48017467491464e-12 2.65328942481270e-12 -5.79570472068573e-13 6.17511954582720e-12 -4.72437536818714e-13
-1.30449341090077e-11 -9.01780276368023e-09 7.92999125348644e-08 4.00936687032002e-11 2.04283104742536e-11 1.24820735380079e-12 9.30760973993229e-12 -7.50079609806551e-12 1.22865516464498e-12 8.06254006030017e-12 5.56252493353000e-12
-1.67007458454893e-12 -4.19265297003500e-11 4.00936687032002e-11 7.92994152491339e-08 -2.20189255359344e-12 -1.19756353847227e-11 9.16324192448759e-12 -8.26074802470566e-12 -2.45588226225298e-12 4.46842214440270e-13 3.39106218157172e-12
-1.69491338826731e-12 -1.21828612766692e-11 2.04283104742536e-11 -2.20189255359344e-12 7.92999951227102e-08 1.69566385897829e-08 -1.13259218230788e-11 -1.06306934205658e-11 5.69865748002221e-13 -1.33330158047386e-11 -2.77602627891578e-12
2.04671335711028e-12 -2.89846125284266e-12 1.24820735380079e-12 -1.19756353847227e-11 1.69566385897829e-08 7.92999831052483e-08 1.16054950164138e-12 5.22206051556241e-12 -7.72954490091433e-12 -2.47645215082267e-12 -3.12393940250136e-12
-5.39085143946466e-12 -2.48017467491464e-12 9.30760973993229e-12 9.16324192448759e-12 -1.13259218230788e-11 1.16054950164138e-12 7.92999555250990e-08 -6.52483339571384e-10 -1.19278810548402e-12 5.25758800878068e-12 6.39317043900566e-12
-6.86316937563370e-12 2.65328942481270e-12 -7.50079609806551e-12 -8.26074802470566e-12 -1.06306934205658e-11 5.22206051556241e-12 -6.52483339571384e-10 7.92994823357396e-08 3.66837619211617e-12 -2.15281354934899e-11 -2.06399667934130e-13
1.13976593936898e-12 -5.79570472068573e-13 1.22865516464498e-12 -2.45588226225298e-12 5.69865748002221e-13 -7.72954490091433e-12 -1.19278810548402e-12 3.66837619211617e-12 7.92997780363045e-08 -3.20478613318390e-09 -9.20728674348331e-13
8.73221272415655e-12 6.17511954582720e-12 8.06254006030017e-12 4.46842214440270e-13 -1.33330158047386e-11 -2.47645215082267e-12 5.25758800878068e-12 -2.15281354934899e-11 -3.20478613318390e-09 7.92999432448874e-08 4.27648531379977e-12
7.29148797499054e-12 -4.72437536818714e-13 5.56252493353000e-12 3.39106218157172e-12 -2.77602627891578e-12 -3.12393940250136e-12 6.39317043900566e-12 -2.06399667934130e-13 -9.20728674348331e-13 4.27648531379977e-12 7.92999522395012e-08];
%% Damping matrix
C=[-3.42898926854769e-06 -3.09853079838774e-10 5.63680179774281e-10 7.18316123505485e-11 7.32043974513657e-11 -8.83641052359976e-11 2.32608057501850e-10 2.96893262006781e-10 -4.94834376380811e-11 -3.77414155634669e-10 -3.15280315155373e-10
-3.09853079838774e-10 -3.42822247644157e-06 3.89849126448742e-07 1.81237812695370e-09 5.26664903076789e-10 1.25368640746015e-10 1.07176084927423e-10 -1.15093176288670e-10 2.53119614277008e-11 -2.66878262734472e-10 2.05557932408960e-11
5.63680179774281e-10 3.89849126448742e-07 -3.42821883018528e-06 -1.73315988046962e-09 -8.82865166143735e-10 -5.37877827811127e-11 -4.02403300820132e-10 3.30719950046673e-10 -5.67461808821010e-11 -3.48413594745865e-10 -2.40423082615653e-10
7.18316123505485e-11 1.81237812695370e-09 -1.73315988046962e-09 -3.42735817917801e-06 9.52596198600402e-11 5.17542525850661e-10 -3.95943650551844e-10 3.56956034992411e-10 1.06211785588847e-10 -1.94110737204741e-11 -1.46516691611563e-10
7.32043974513657e-11 5.26664903076789e-10 -8.82865166143735e-10 9.52596198600402e-11 -3.42682637747221e-06 -7.32754904764228e-07 4.89474123375667e-10 4.59411557481755e-10 -2.46702162282727e-11 5.75962888883721e-10 1.19725151113121e-10
-8.83641052359976e-11 1.25368640746015e-10 -5.37877827811127e-11 5.17542525850661e-10 -7.32754904764228e-07 -3.42682555369805e-06 -5.00861421182035e-11 -2.25669451568986e-10 3.33860804725077e-10 1.07182548547866e-10 1.34689134354363e-10
2.32608057501850e-10 1.07176084927423e-10 -4.02403300820132e-10 -3.95943650551844e-10 4.89474123375667e-10 -5.00861421182035e-11 -3.42599149801318e-06 2.81892081456211e-08 5.16487620733699e-11 -2.26966026963761e-10 -2.76288740535447e-10
2.96893262006781e-10 -1.15093176288670e-10 3.30719950046673e-10 3.56956034992411e-10 4.59411557481755e-10 -2.25669451568986e-10 2.81892081456211e-08 -3.42597017380682e-06 -1.50581477924180e-10 9.29819231192795e-10 9.05872401930873e-12
-4.94834376380811e-11 2.53119614277008e-11 -5.67461808821010e-11 1.06211785588847e-10 -2.46702162282727e-11 3.33860804725077e-10 5.16487620733699e-11 -1.50581477924180e-10 -3.42487796344742e-06 1.38411425482107e-07 3.97409661925509e-11
-3.77414155634669e-10 -2.66878262734472e-10 -3.48413594745865e-10 -1.94110737204741e-11 5.75962888883721e-10 1.07182548547866e-10 -2.26966026963761e-10 9.29819231192795e-10 1.38411425482107e-07 -3.42488187253882e-06 -1.84679726986698e-10
-3.15280315155373e-10 2.05557932408960e-11 -2.40423082615653e-10 -1.46516691611563e-10 1.19725151113121e-10 1.34689134354363e-10 -2.76288740535447e-10 9.05872401930873e-12 3.97409661925509e-11 -1.84679726986698e-10 -3.42460989035060e-06];
%% Stiffness matrix
K=[0.763518187446834 -0.123471983781582 2.12644497110045 -0.00159184932062271 0.00286665123297165 0.00326795662599221 -0.00565823204629348 -0.0111247074471964 -0.852444758506581 0.0170129698870622 0.000567993613223703
0.122986250969018 4.76907928529478 -0.543239600079982 2.12608568267782 0.220268589553943 0.0633458824464058 -0.00720500630476635 5.67455459314993e-05 0.000218541703671969 0.00504540957654056 0.00297905359335030
-2.13057558844554 -0.544693846047297 4.79735395984435 -0.116535843976208 -3.80931536422324 -1.08160690561449 0.00135336667261543 0.00474788044688154 0.00305299632669883 0.00159367036253352 -0.0179988437727940
-0.00372848570465477 -2.13382763678292 0.123763968250996 12.3142831297151 -0.000137283537755157 -0.00228005296293815 3.81673953925391 -0.0591314258833513 0.00127018393509301 -0.00796895715491763 0.00234941752038280
-0.00163737934937908 -0.222001483568310 3.81734112223475 0.00173003724327611 19.0472196037760 4.07106712942156 -0.00255853650434496 -0.125116265764028 5.46695676693670 -0.0993709375445011 0.00703634900851892
0.00258298057052472 -0.0623605607181842 1.08529817957274 -0.00267072015732558 4.07294423066431 19.0574340781544 0.0147753791569810 2.05513925508344 1.52210287227782 -0.0237523853866786 0.00922809527607430
-0.00468535654428222 -0.00626853319266933 -0.00115469998983706 -3.81119446621267 -0.00155959499500691 -0.0121660190928631 32.1310954783617 -0.262056094807989 0.000541439656752278 0.00464738145837403 -0.0255209825354020
0.0134976568920513 -0.000376978112019705 0.00263963824092268 0.0542232932656596 0.120732014828162 -2.05279502259563 -0.263280871097288 32.3912983853784 0.00480792646303815 -0.0114157922244635 -3.80017279650533
0.847260618562458 0.000667255021819536 0.00160008219803174 0.000312861305428450 -5.46833245082724 -1.53174909321691 0.00133283528195683 0.00562883566011632 55.0224400988869 -2.22476909528779 0.00148357487599929
-0.0121173666954115 0.00477027008481583 0.00127363569963892 0.00477817239843339 0.0854464228914320 0.0287153890552796 0.00377260121504029 -0.0113582049525324 -2.22603066376352 55.0972049682465 0.00355018925501934
0.000736959774981387 0.00110272425151785 0.0202945683248242 0.000494668031994278 0.00591913812055584 0.00913408539493314 0.0270935852755023 3.80595460544076 -0.00385777038321837 0.00320054596414318 61.4845746182663];
%% Force vector
F=[-0.00099996606222035952475160617808552
-0.0000002460836756585278896861534876451
-0.00000060555830113300843168099357979169
-0.00000023236067999397214860565123562071
-0.00092724784232042312137180967502376
-0.0012336739219426489876333403630418
0.00000015739906629771823481279895519369
0.000000049705400083053920707370130069998
0.0000000074445232833242759638623669833003
0.0000001051612538093116976310734137023
-0.00099218997013689166348230850119706];
%% State matrix
A = [zeros(size(M)), diag(ones(1, length(M))); -M\K, -M\C];
%% Input matrix
B = [zeros(size(M)); diag(ones(1, length(M)))];
%% State equations (Matrix form)
dvdt = A*v + B*F;
end
More Answers (2)
syms t tau_1(t) tau_2(t) tau_3(t) tau_4(t) tau_5(t) tau_6(t) tau_7(t) tau_8(t) tau_9(t) tau_10(t) tau_11(t)
v = transpose([tau_1(t) tau_2(t) tau_3(t) tau_4(t) tau_5(t) tau_6(t) tau_7(t) tau_8(t) tau_9(t) tau_10(t) tau_11(t)])
odes = M * diff(diff(v,t),t) + C * diff(v,t) + K * v - F == 0;
[V,S] = odeToVectorField(odes)
fun_ode = matlabFunction(V,'vars',{'t','Y'});
Y0 = zeros(2*numel(v),1);
tspan = [0, 1];
[t, y] = ode15s(fun_ode, tspan, Y0);
plot(t, y)
Hassaan
on 4 Jan 2024
% Define matrices M, C, K, and vector F (Use your actual matrices and vector here)
M = [...];
C = [...];
K = [...];
F = [...];
% Define the symbolic variables
syms tau_1(t) tau_2(t) tau_3(t) tau_4(t) tau_5(t) tau_6(t) tau_7(t) tau_8(t) tau_9(t) tau_10(t) tau_11(t)
% Define the state vector and its derivative
v = [tau_1; tau_2; tau_3; tau_4; tau_5; tau_6; tau_7; tau_8; tau_9; tau_10; tau_11];
dv = diff(v, t);
ddv = diff(dv, t);
% Define the ODEs
odes = M * ddv + C * dv + K * v - F;
% Convert the symbolic ODEs to a MATLAB function handle
V = matlabFunction(odes, 'Vars', {t, [v; dv]});
% Define initial conditions for v and dv
Y0 = zeros(22, 1); % Replace with your actual initial conditions
% Define the time span
tspan = [0, 1];
% Create a wrapper function for ode45
odefun = @(t, Y) [Y(12:end); V(t, Y)];
% Solve the ODE system
[t, y] = ode45(odefun, tspan, Y0);
% Plot the results
plot(t, y(:, 1:11)); % Plot only the positions (first 11 state variables)
Remember, when using ode45, you need to convert the second-order ODEs into a set of first-order ODEs by introducing new variables to represent the velocities, which it seems you have tried to do with dv. The plot command at the end is tailored to show the results of the first 11 state variables, which represent the tau_i variables you have in your system.
------------------------------------------------------------------------------------------------------------------------------------------------
If you find the solution helpful and it resolves your issue, it would be greatly appreciated if you could accept the answer. Also, leaving an upvote and a comment are also wonderful ways to provide feedback.
Professional Interests
- Technical Services and Consulting
- Embedded Systems | Firmware Development | Simulations
- Electrical and Electronics Engineering
Categories
Find more on Mathematics in Help Center and File Exchange
