solving a system of 2nd order 2 ODE: output vector length not matching the initial condition vector length

I am trying to solve a 2nd order system of 2 ODEs using the ODE45. the variables to be solved are V_1 and V_6.

the output vetor y = [V_1 ; V_1' ; V_6 ;V_6'] as I defined has the length of 4, and it matches the definition of the system of ODEs

the length of the initial condition vector is 4, and it matches the definition of the system of the 2 ODEs

but I get this error:

"LINEAR_ODE returns a vector of length 2, but the length of initial conditions vector is 4. The vector returned by

LINEAR_ODE and the initial conditions vector must have the same number of elements."

the ODE equations are pretty complicated but correspond to the equations of a coupled two LC resonators. I checked my code pretty much but going crazy.

I am following the same approch as is in this question (https://www.mathworks.com/matlabcentral/answers/491409-how-to-solve-system-of-2nd-order-differential-equations-using-ode45), and the code works with no error !

Thank you so much !

TimeScale= 10*1E-10;     
TimeStep= 1 * 1E-15;

tspan = [0:TimeStep: TimeScale];

y0 = [0; 0; 0;0];

% y = [V_1; V_1'; V_6 ;V_6'] 
%y(1)=V_1  y(2)=V_1'  y(3)=V_6 y(4)=V_6'

[t, y]= ode45(@Linear_ODE, tspan, y0);


function dydt=Linear_ODE(t,y)

h = 6.62607015 * 1E-34 ;  % [J s]
e=1.60217662 * 1E-19;     %[C]   % elementary charge(C)
K_J= (2*e*2*pi)/h ;%2e/h_bar

C_1 = 0.3e-12;

L_5 = 10.0e-9;

L_J = 2*L_5;
L_3 = L_J; 

M_4 = 100e-12;

C_6 = 0.3e-12;

tildeL_3 = L_3-M_4; 
tildeL_5 = L_5-M_4;

V_0 = 1e-6;
R_s=50;

Ic_0= 25e-9;

w= 2.9e09;
Alpha= tildeL_3 + M_4;
Beta = tildeL_5 + M_4; 

TimeStep= 1 * 1E-15;

dydt= [y(2); -1/(R_s*C_1)*y(2)-(1 + Alpha * 1/L_J)*1/(Alpha*C_1) * y(1) + M_4*C_6/(C_1*Alpha) * diff(y(4))/TimeStep-V_0*w*sin(w*t)/(R_s*C_1) ; y(4) ; M_4*C_1/(C_6*Beta)* diff(y(2))/TimeStep+ M_4/(R_s*C_6*Beta)* y(2)+ M_4*1/(C_6*Beta)* 1/L_J * y(1)+y(3)/(C_6*Beta)+ M_4/(R_s*C_6*Beta)*V_0*w*cos(w*t)];


end

Answer 1

[David Goodmanson]

Hi yosri,

the reason this is happening is that the expression for dy/dt has four elements, but one of them contains diff(y(4)) and another contains diff(y(2)). Since y(4) is a scalar, diff(y(4)) = [ ], i.e. null. Same for y(2). The nulls collapse the dydt vector from four elements down to 2, hence the error.

Right now you appear to have a second order differential equation in two variables (call them a and b), so that

y(1) = a   y(2) = a'  y(3) = b   y(4) = b' 

It is not clear what the intent of diff(y(4)) is, but ordinarily you would be writing out dydt(4) = b'' as a function of time and the y variables, with no diffs involved.

I guess it's personal preference, but long lines of code like the one for dydt do not lend themselves well to error checking. It doesn't hurt to at least break up the line into sections with the use of dotdotdot.

[yosri ayadi]

Hi David,

Thank you so much! I understand the mistake now.

I have two coupled ODEs in this form let's say:

a" = alpha a' + Beta a + Gamma b"

b"= theta b'+ epsilon b+ psi a"

In order to represent the b" term I tried: diff(y(3),t,2) and then diff(y(4))/TimeStep

same thing to represent a"

I thought that y(4) is a vector along with the integration...

Do you advise on anyway to represent the a" term ?

thx a lot !

[David Goodmanson]

Hi yosri

Consider a,a',b,b' to be known. You have two linear equations in the two unknowns, a'' and b''. Solve those two equations to find each of a'' and b'' in terms of a,a',b,b', and then go back to ode45.

[yosri ayadi]

do you mean solving using the symbolic sys ?

I am solving the two coupled ODEs to find a and b. the issue is that in both ODEs a" and b" are involved in the equation, i.e. a"=fct(a, a', b") and b"=(b, b', a").

thank you

[David Goodmanson]

if you consider everything except a'' and b'' to be constants, then you have two linear equations in the two unknowns a'' and b''. You just solve those in linear algebra fashion to find a'' and b'' as a function of the 'constants'. You should be able to do that by hand. Or, if for notational purposes you get rid of the quote marks by renaming a' as a1, a'' as a2 etc (or whatever), then you can get the symbolic solver to do it.