Problem in solving matrix diiferential equation using ode45?

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 ;

Error using indexing
Invalid indexing or function definition. Indexing must follow MATLAB indexing. Function arguments must be symbolic variables, and function body must be sym expression.

Error in indexing (line 968)
R_tilde = builtin('subsref',L_tilde,Idx);

Y = matlabFunction(odes,'Vars',v(t));
Y0     = [0;0];
tspan  = [0, 1];
[t, y] = ode45(@Y, tspan, Y0);
plot(t, y);

[Dyuman Joshi]

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.

You need to use odeToVectorField and simplify the equations, then proceed.

Edit - See @Torsten's answer.

Answer 1

[Sam Chak]

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

Answer 2

[Torsten]

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)

Answer 3

[Hassaan]

% 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