Trying to use Symbolic Math Toolbox to prepare to solve two coupled second order ODEs

Hello, I have the following 2nd order differential equations I need to solve:

Fex1 = (m1+A11)*z1''+A12*z2''+(b11+bvis1)*z1'+b12*z2'+(b+bmech)*(z1'−z2')+C*(z1−z2)+Cb*z1

Fex2 = (m2+A22)*z2''+A21*z1''+(b22+bvis2)*z2'+b21*z1'+(b+bmech)*(z2'−z1')+C*(z2−z1)

When googling about how to go about it i stumbled upon a Post from 2 years ago that seemed to have a pretty similar problem to solve. https://www.mathworks.com/matlabcentral/answers/398099-use-ode45-to-solve-a-system-of-two-coupled-second-order-odes

my Parameters, not all of them are known yet, so most of them are random values to simply get the code to working.

%Parameter
r = 9.5; 
rho_w = 0.001; 
g = 9.81; 

m_1 = 0.913; 
m_2p = 0.79;
m_2e = 1.2; 
m_2 = m_2p+m_2e; 
A_11 = 0.1; 
A_12 = 0.1; 
A_21 = 0.1; 
A_22 = 0.1; 
b = 0.1; 
b_mech = 0.002806;
b_vis1 = 0.1; 
b_vis2 = 0.1; 
b_11 = 0.1; 
b_12 = 0.1; 
b_21 = 0.1; 
b_22 = 0.1; 
C = 0.25; 
Cb = rho_w*g*pi*r^2; 

A = 0.1; 
w = 3; 
tspan=[0:0.01:25]; 
X_1 = 15;
X_2 = 5;
Ftfcn1 = @(t) A*X_1*sin(w*t);
Ftfcn2 = @(t) A*X_2*sin(w*t);

So i tried to solve it the way they did.

syms m_1 m_2 A_11 A_12 A_21 A_22 b b_11 b_12 b_21 b_22 b_vis1 b_vis2 b_mech C Cb Ftfcn1 Ftfcn2 z1(t) z2(t) Y
%z1'' =
%(F_ex1-A_12*z2''-(b_11+b_vis1)*z1'-b_12*z2'-(b+b_mech)*(z1'-z2')-C*(z1-z2)-Cb*z1)/(m_1+A_11)
Dz1 = diff(z1);
D2z1 = diff(z1,2);
Dz2 = diff(z2);
D2z2 = diff(z2,2);
Eq1 = D2z1 == (Ftfcn1-A_12*D2z2-(b_11+b_vis1)*Dz1-b_12*Dz2-(b+b_mech)*(Dz1-Dz2)-C*(z1-z2)-Cb*z1)/(m_1+A_11);
%z2'' =
%(F_ex2-A_21*z1''-(b_22+b_vis2)*z2'-b_21*z1'-(b+b_mech)*(z2'-z1')-C*(z2-z1))/(m_2+A_22)
Eq2 = D2z2 == (Ftfcn2-A_21*D2z1-(b_22+b_vis2)*Dz2-b_21*Dz1-(b+b_mech)*(Dz2-Dz1)-C*(z2-z1))/(m_2+A_22);
[VF,Subs] = odeToVectorField(Eq1,Eq2);
ftotal = matlabFunction(VF, 'Vars',{t,Y,m_1, m_2, A_11, A_12, A_21, A_22, b, b_11, b_12, b_21, b_22, b_vis1, b_vis2, b_mech, C, Cb, Ftfcn1, Ftfcn2});

After edditing the ftotal to include the (t) for the Force functions i ended up with (quite long)

ftotal = @(t,Y,m_1,m_2,A_11,A_12,A_21,A_22,b,b_11,b_12,b_21,b_22,b_vis1,b_vis2,b_mech,C,Cb,Ftfcn1,Ftfcn2)[Y(2);-(-Ftfcn2(t).*m_1-A_11.*Ftfcn2(t)+A_21.*Ftfcn1(t)+A_11.*C.*Y(1)-A_11.*C.*Y(3)+A_21.*C.*Y(1)-A_21.*C.*Y(3)-A_21.*Cb.*Y(3)+A_11.*b.*Y(2)-A_11.*b.*Y(4)+A_21.*b.*Y(2)-A_21.*b.*Y(4)+A_11.*b_22.*Y(2)-A_21.*b_12.*Y(2)+A_11.*b_21.*Y(4)-A_21.*b_11.*Y(4)+A_11.*b_mech.*Y(2)-A_11.*b_mech.*Y(4)+A_21.*b_mech.*Y(2)-A_21.*b_mech.*Y(4)+A_11.*b_vis2.*Y(2)-A_21.*b_vis1.*Y(4)+C.*m_1.*Y(1)-C.*m_1.*Y(3)+b.*m_1.*Y(2)-b.*m_1.*Y(4)+b_22.*m_1.*Y(2)+b_21.*m_1.*Y(4)+b_mech.*m_1.*Y(2)-b_mech.*m_1.*Y(4)+b_vis2.*m_1.*Y(2))./(A_11.*m_2+A_22.*m_1+m_1.*m_2+A_11.*A_22-A_12.*A_21);Y(4);-(-Ftfcn1(t).*m_2+A_12.*Ftfcn2(t)-A_22.*Ftfcn1(t)-A_12.*C.*Y(1)+A_12.*C.*Y(3)-A_22.*C.*Y(1)+A_22.*C.*Y(3)+A_22.*Cb.*Y(3)-A_12.*b.*Y(2)+A_12.*b.*Y(4)-A_22.*b.*Y(2)+A_22.*b.*Y(4)-A_12.*b_22.*Y(2)+A_22.*b_12.*Y(2)-A_12.*b_21.*Y(4)+A_22.*b_11.*Y(4)-A_12.*b_mech.*Y(2)+A_12.*b_mech.*Y(4)-A_22.*b_mech.*Y(2)+A_22.*b_mech.*Y(4)-A_12.*b_vis2.*Y(2)+A_22.*b_vis1.*Y(4)-C.*m_2.*Y(1)+C.*m_2.*Y(3)+Cb.*m_2.*Y(3)-b.*m_2.*Y(2)+b.*m_2.*Y(4)+b_12.*m_2.*Y(2)+b_11.*m_2.*Y(4)-b_mech.*m_2.*Y(2)+b_mech.*m_2.*Y(4)+b_vis1.*m_2.*Y(4))./(A_11.*m_2+A_22.*m_1+m_1.*m_2+A_11.*A_22-A_12.*A_21)]

So for my ode45 function

ic = zeros(4,1);
ic(1,1) = A;
[T,Y] = ode45(@(t,Y) ftotal(t,Y,m_1,m_2,A_11,A_12,A_21,A_22,b,b_11,b_12,b_21,b_22,b_vis1,b_vis2,b_mech,C,Cb,Ftfcn1,Ftfcn2), tspan, ic);
plot(T, Y)
grid

I have obviously mucked something up since im getting quite a few error codes (not sure if they help or not, but here they are):

Array indices must be positive integers or logical values.

Error in sym/subsref (line 902)
            R_tilde = builtin('subsref',L_tilde,Idx);

Error in
Bewegungsgleichung_mit_Ode45>@(t,Y,m_1,m_2,A_11,A_12,A_21,A_22,b,b_11,b_12,b_21,b_22,b_vis1,b_vis2,b_mech,C,Cb,Ftfcn1,Ftfcn2)[Y(2);-(-Ftfcn2(t).*m_1-A_11.*Ftfcn2(t)+A_21.*Ftfcn1(t)+A_11.*C.*Y(1)-A_11.*C.*Y(3)+A_21.*C.*Y(1)-A_21.*C.*Y(3)-A_21.*Cb.*Y(3)+A_11.*b.*Y(2)-A_11.*b.*Y(4)+A_21.*b.*Y(2)-A_21.*b.*Y(4)+A_11.*b_22.*Y(2)-A_21.*b_12.*Y(2)+A_11.*b_21.*Y(4)-A_21.*b_11.*Y(4)+A_11.*b_mech.*Y(2)-A_11.*b_mech.*Y(4)+A_21.*b_mech.*Y(2)-A_21.*b_mech.*Y(4)+A_11.*b_vis2.*Y(2)-A_21.*b_vis1.*Y(4)+C.*m_1.*Y(1)-C.*m_1.*Y(3)+b.*m_1.*Y(2)-b.*m_1.*Y(4)+b_22.*m_1.*Y(2)+b_21.*m_1.*Y(4)+b_mech.*m_1.*Y(2)-b_mech.*m_1.*Y(4)+b_vis2.*m_1.*Y(2))./(A_11.*m_2+A_22.*m_1+m_1.*m_2+A_11.*A_22-A_12.*A_21);Y(4);-(-Ftfcn1(t).*m_2+A_12.*Ftfcn2(t)-A_22.*Ftfcn1(t)-A_12.*C.*Y(1)+A_12.*C.*Y(3)-A_22.*C.*Y(1)+A_22.*C.*Y(3)+A_22.*Cb.*Y(3)-A_12.*b.*Y(2)+A_12.*b.*Y(4)-A_22.*b.*Y(2)+A_22.*b.*Y(4)-A_12.*b_22.*Y(2)+A_22.*b_12.*Y(2)-A_12.*b_21.*Y(4)+A_22.*b_11.*Y(4)-A_12.*b_mech.*Y(2)+A_12.*b_mech.*Y(4)-A_22.*b_mech.*Y(2)+A_22.*b_mech.*Y(4)-A_12.*b_vis2.*Y(2)+A_22.*b_vis1.*Y(4)-C.*m_2.*Y(1)+C.*m_2.*Y(3)+Cb.*m_2.*Y(3)-b.*m_2.*Y(2)+b.*m_2.*Y(4)+b_12.*m_2.*Y(2)+b_11.*m_2.*Y(4)-b_mech.*m_2.*Y(2)+b_mech.*m_2.*Y(4)+b_vis1.*m_2.*Y(4))./(A_11.*m_2+A_22.*m_1+m_1.*m_2+A_11.*A_22-A_12.*A_21)]
(line 107)
ftotal =
@(t,Y,m_1,m_2,A_11,A_12,A_21,A_22,b,b_11,b_12,b_21,b_22,b_vis1,b_vis2,b_mech,C,Cb,Ftfcn1,Ftfcn2)[Y(2);-(-Ftfcn2(t).*m_1-A_11.*Ftfcn2(t)+A_21.*Ftfcn1(t)+A_11.*C.*Y(1)-A_11.*C.*Y(3)+A_21.*C.*Y(1)-A_21.*C.*Y(3)-A_21.*Cb.*Y(3)+A_11.*b.*Y(2)-A_11.*b.*Y(4)+A_21.*b.*Y(2)-A_21.*b.*Y(4)+A_11.*b_22.*Y(2)-A_21.*b_12.*Y(2)+A_11.*b_21.*Y(4)-A_21.*b_11.*Y(4)+A_11.*b_mech.*Y(2)-A_11.*b_mech.*Y(4)+A_21.*b_mech.*Y(2)-A_21.*b_mech.*Y(4)+A_11.*b_vis2.*Y(2)-A_21.*b_vis1.*Y(4)+C.*m_1.*Y(1)-C.*m_1.*Y(3)+b.*m_1.*Y(2)-b.*m_1.*Y(4)+b_22.*m_1.*Y(2)+b_21.*m_1.*Y(4)+b_mech.*m_1.*Y(2)-b_mech.*m_1.*Y(4)+b_vis2.*m_1.*Y(2))./(A_11.*m_2+A_22.*m_1+m_1.*m_2+A_11.*A_22-A_12.*A_21);Y(4);-(-Ftfcn1(t).*m_2+A_12.*Ftfcn2(t)-A_22.*Ftfcn1(t)-A_12.*C.*Y(1)+A_12.*C.*Y(3)-A_22.*C.*Y(1)+A_22.*C.*Y(3)+A_22.*Cb.*Y(3)-A_12.*b.*Y(2)+A_12.*b.*Y(4)-A_22.*b.*Y(2)+A_22.*b.*Y(4)-A_12.*b_22.*Y(2)+A_22.*b_12.*Y(2)-A_12.*b_21.*Y(4)+A_22.*b_11.*Y(4)-A_12.*b_mech.*Y(2)+A_12.*b_mech.*Y(4)-A_22.*b_mech.*Y(2)+A_22.*b_mech.*Y(4)-A_12.*b_vis2.*Y(2)+A_22.*b_vis1.*Y(4)-C.*m_2.*Y(1)+C.*m_2.*Y(3)+Cb.*m_2.*Y(3)-b.*m_2.*Y(2)+b.*m_2.*Y(4)+b_12.*m_2.*Y(2)+b_11.*m_2.*Y(4)-b_mech.*m_2.*Y(2)+b_mech.*m_2.*Y(4)+b_vis1.*m_2.*Y(4))./(A_11.*m_2+A_22.*m_1+m_1.*m_2+A_11.*A_22-A_12.*A_21)]

Error in Bewegungsgleichung_mit_Ode45>@(t,Y)ftotal(t,Y,m_1,m_2,A_11,A_12,A_21,A_22,b,b_11,b_12,b_21,b_22,b_vis1,b_vis2,b_mech,C,Cb,Ftfcn1,Ftfcn2) (line 112)
[T,Y] = ode45(@(t,Y) ftotal(t,Y,m_1,m_2,A_11,A_12,A_21,A_22,b,b_11,b_12,b_21,b_22,b_vis1,b_vis2,b_mech,C,Cb,Ftfcn1,Ftfcn2), tspan, ic);

Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)

I am totally out of my comfortzone here. I have tried to get this to working for a while, but i am unable to do so. I dont think my initial functions are incorrect, but i wouldnt bet my life on it. Is the Symbolic Math Toolbox a good way to go about it ? It looked to easy and the code in the old question works perfectly. Thank you if anyone can help me here.

[Star Strider]

What are ‘Fex1’ and ‘Fex2’ derivatives of?

[Benedict Bluhm]

Sorry i should have really not asked this at night right before going to sleep.

Fex1 = (m1+A11)*z1''+A12*z2''+(b11+bvis1)*z1'+b12*z2'+(b+bmech)*(z1'−z2')+C*(z1−z2)+Cb*z1

Fex2 = (m2+A22)*z2''+A21*z1''+(b22+bvis2)*z2'+b21*z1'+(b+bmech)*(z2'−z1')+C*(z2−z1)

should be rearranged to

z1'' = (Fex1(t) - A12*z2''-(b11+bvis1)*z1'-b12*z2'-(b+bmech)*(z1'−z2')-C*(z1−z2)-Cb*z1)/(m1+A11)

z2'' = (Fex2(t) - A21*z1''-(b22+bvis2)*z2'-b21*z1'-(b+bmech)*(z2'−z1')-C*(z2−z1))/(m2+A22)

Derivatives are with respect to time.

I renamed Fex1/2 to Ftfcn1/2 in the code.

I just realised that you are the one that originally replied to the other Question two years ago. Sorry if i massacred the code you wrote.

Answer 1

[Walter Roberson]

After edditing the ftotal to include the (t) for the Force functions

Not a good idea.

A = 0.1; 
w = 3; 
X_1 = 15;
X_2 = 5;
Ftfcn1 = @(t) A*X_1*sin(w*t);
Ftfcn2 = @(t) A*X_2*sin(w*t);

syms m_1 m_2 A_11 A_12 A_21 A_22 b b_11 b_12 b_21 b_22 b_vis1 b_vis2 b_mech C Cb z1(t) z2(t) Y
Dz1 = diff(z1);
D2z1 = diff(z1,2);
Dz2 = diff(z2);
D2z2 = diff(z2,2);
Eq1 = D2z1 == (Ftfcn1(t)-A_12*D2z2-(b_11+b_vis1)*Dz1-b_12*Dz2-(b+b_mech)*(Dz1-Dz2)-C*(z1-z2)-Cb*z1)/(m_1+A_11);
Eq2 = D2z2 == (Ftfcn2(t)-A_21*D2z1-(b_22+b_vis2)*Dz2-b_21*Dz1-(b+b_mech)*(Dz2-Dz1)-C*(z2-z1))/(m_2+A_22);
[VF,Subs] = odeToVectorField(Eq1,Eq2);
ftotal = matlabFunction(VF, 'Vars',{t,Y,m_1, m_2, A_11, A_12, A_21, A_22, b, b_11, b_12, b_21, b_22, b_vis1, b_vis2, b_mech, C, Cb})

ftotal = function_handle with value:
    @(t,Y,m_1,m_2,A_11,A_12,A_21,A_22,b,b_11,b_12,b_21,b_22,b_vis1,b_vis2,b_mech,C,Cb)[Y(2);((-m_1.*sin(t.*3.0)-A_11.*sin(t.*3.0)+A_21.*sin(t.*3.0).*3.0+A_11.*C.*Y(1).*2.0-A_11.*C.*Y(3).*2.0+A_21.*C.*Y(1).*2.0-A_21.*C.*Y(3).*2.0-A_21.*Cb.*Y(3).*2.0+A_11.*b.*Y(2).*2.0-A_11.*b.*Y(4).*2.0+A_21.*b.*Y(2).*2.0-A_21.*b.*Y(4).*2.0+A_11.*b_22.*Y(2).*2.0-A_21.*b_12.*Y(2).*2.0+A_11.*b_21.*Y(4).*2.0-A_21.*b_11.*Y(4).*2.0+A_11.*b_mech.*Y(2).*2.0-A_11.*b_mech.*Y(4).*2.0+A_21.*b_mech.*Y(2).*2.0-A_21.*b_mech.*Y(4).*2.0+A_11.*b_vis2.*Y(2).*2.0-A_21.*b_vis1.*Y(4).*2.0+C.*m_1.*Y(1).*2.0-C.*m_1.*Y(3).*2.0+b.*m_1.*Y(2).*2.0-b.*m_1.*Y(4).*2.0+b_22.*m_1.*Y(2).*2.0+b_21.*m_1.*Y(4).*2.0+b_mech.*m_1.*Y(2).*2.0-b_mech.*m_1.*Y(4).*2.0+b_vis2.*m_1.*Y(2).*2.0).*(-1.0./2.0))./(A_11.*m_2+A_22.*m_1+m_1.*m_2+A_11.*A_22-A_12.*A_21);Y(4);((m_2.*sin(t.*3.0).*-3.0+A_12.*sin(t.*3.0)-A_22.*sin(t.*3.0).*3.0-A_12.*C.*Y(1).*2.0+A_12.*C.*Y(3).*2.0-A_22.*C.*Y(1).*2.0+A_22.*C.*Y(3).*2.0+A_22.*Cb.*Y(3).*2.0-A_12.*b.*Y(2).*2.0+A_12.*b.*Y(4).*2.0-A_22.*b.*Y(2).*2.0+A_22.*b.*Y(4).*2.0-A_12.*b_22.*Y(2).*2.0+A_22.*b_12.*Y(2).*2.0-A_12.*b_21.*Y(4).*2.0+A_22.*b_11.*Y(4).*2.0-A_12.*b_mech.*Y(2).*2.0+A_12.*b_mech.*Y(4).*2.0-A_22.*b_mech.*Y(2).*2.0+A_22.*b_mech.*Y(4).*2.0-A_12.*b_vis2.*Y(2).*2.0+A_22.*b_vis1.*Y(4).*2.0-C.*m_2.*Y(1).*2.0+C.*m_2.*Y(3).*2.0+Cb.*m_2.*Y(3).*2.0-b.*m_2.*Y(2).*2.0+b.*m_2.*Y(4).*2.0+b_12.*m_2.*Y(2).*2.0+b_11.*m_2.*Y(4).*2.0-b_mech.*m_2.*Y(2).*2.0+b_mech.*m_2.*Y(4).*2.0+b_vis1.*m_2.*Y(4).*2.0).*(-1.0./2.0))./(A_11.*m_2+A_22.*m_1+m_1.*m_2+A_11.*A_22-A_12.*A_21)]

[Benedict Bluhm]

Could you explain why? Im fairly new to this and it worked when i trie it for the different code i linked, which is why tried to replicate it for myself.

[Walter Roberson]

We can be relatively sure that the problem you had is that after you edited ftotal to include the (t), that you failed to do

Ftfcn1 = @(t) A*X_1*sin(w*t);
Ftfcn2 = @(t) A*X_2*sin(w*t);

before the call to

[T,Y] = ode45(@(t,Y) ftotal(t,Y,m_1,m_2,A_11,A_12,A_21,A_22,b,b_11,b_12,b_21,b_22,b_vis1,b_vis2,b_mech,C,Cb,Ftfcn1,Ftfcn2), tspan, ic);

and as a result were passing in the still-symbolic scalars Ftfcn1 and Ftfcn2 instead of the force functions.

Hand editing the results of matlabFunction is prone to error.

[Benedict Bluhm]

Thank you that makes sense. How would you solve the equations without hand editing the ftotal ?

[Walter Roberson]

The way I already posted at https://www.mathworks.com/matlabcentral/answers/719490-trying-to-use-symbolic-math-toolbox-to-prepare-to-solve-two-coupled-second-order-odes#answer_600000 does not require hand-editing the ftotal.

[Benedict Bluhm]

Thank you very much. That has cut down the amount of errors i get by quite a bit.

Error using odearguments (line 113)
Inputs must be floats, namely single or double.

Error in ode45 (line 115)
  odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in BewegungsgleichungToolBox (line 52)
[T,Y] = ode45(@(t,Y) ftotal(t,Y,m_1,m_2,A_11,A_12,A_21,A_22,b,b_11,b_12,b_21,b_22,b_vis1,b_vis2,b_mech,C,Cb), tspan, ic);

Obviously something isnt quite right yet, but im happy about any progress currently.

[Walter Roberson]

Either your tspan or your ic contains symbolic values, possibly.

The ode*() routines are strict that the tspan and initial conditions and results of the function are all single or double class, and will complain about anything else, even if the other thing is convertable to double.

[Benedict Bluhm]

Thank you so much i finally figured it out. Youve been a big help.