How do I solve a system of differential equations with initial conditions to the second derivative?

Hi! I'm kind of a newby regarding Matlab, and I have run into trouble trying to solve a modelling problem given by my university. I really hope I can get some help here since I'm quite obsessed by finding a decent solution.

I have to model a double pendulum with certain initial conditions in Matlab. For this, I set up some equations with symbolic variables in them and differentiated some of these to get a system containing second order differential equations. I also have some initial conditions that need to be applied. These initial conditions regard the initial symbolic variables and their first derivatives, so the unknowns of the functions have now become the second derivatives of the initial symbolic variables. I thus have to solve the system of equations, including the constraints, for these second derivatives. This puzzles me, since I can't seem to find a way to do this.

I have been searching this forum and google in general for almost two hours now, and I'm starting to get desparate! I already found a way to solve systems of differential equations using dsolve ( this link ) as well as how to include initial conditions ( this link ), but I can't seem to find how to tell the solver that I want to solve for the second derivatives of my symbolic variables. I already tried just running it but Matlab goes into total lockdown and never stops computating, so I had to brute-force a restart to make it workable again.

As was my story, the code regarding this problem is quite extensive but I will post the entire thing anyway:

%%Initiation
% Quantities
rho1 = 1180; %kg/m^3

rho2 = rho1;
L1 = 0.55; %m
L2 = L1;
w1 = 50/1000; %m
w2 = w1;
t1 = 4/1000; %m
t2 = t1;
g = 9.81; %N/kg

m1 = w1*t1*L1*rho1;
m2 = w2*t2*L2*rho2;
I1 = 1/12*m1*(w1^2+t1^2);
I2 = 1/12*m2*(w2^2+t2^2);

% Unknowns
syms x1(t) y1(t) theta1(t) x2(t) y2(t) theta2(t) HA VA HB VB

dx1 = diff(x1);
dy1 = diff(y1);
dtheta1 = diff(theta1);
dx2 = diff(x2);
dy2 = diff(y2);
dtheta2 = diff(theta2);

% ddx1 = diff(x1,2);
% ddy1 = diff(y1,2);
% ddtheta1 = diff(theta1,2);
% ddx2 = diff(x2,2);
% ddy2 = diff(y2,2);
% ddtheta2 = diff(theta2,2);

%%a.
% Equations of motion:
% Body 1
H1 = m1*diff(x1,2);
V1 = m1*diff(y1,2);
M1 = I1*diff(theta1,2);

motion_1x = HA - HB + m1*g == H1;
motion_1y = VA - VB == V1;
motion_1m = (HA + HB)*sin(theta1)*L1/2 - (VA + VB)*cos(theta1)*L1/2 == M1;

% Body 2
H2 = m2*diff(x2,2);
V2 = m2*diff(y2,2);
M2 = I2*diff(theta2,2);

motion_2x = HB + m2*g == H2;
motion_2y = VB == V2;
motion_2m = HB*sin(theta2)*L2/2 - VB*cos(theta2)*L2/2 == M2;

% Constraint equations:
% Joint A
constr_Ax = x1 - cos(theta1)*L1/2 == 0;
constr_Ay = y1 - sin(theta1)*L1/2 == 0;

constr_Ax = diff(constr_Ax,2);
constr_Ay = diff(constr_Ay,2);

% Joint B
constr_Bx = x1 + cos(theta1)*L1/2 == x2 - cos(theta2)*L2/2;
constr_By = y1 + sin(theta1)*L1/2 == y2 - sin(theta2)*L2/2;

constr_Bx = diff(constr_Bx,2);
constr_By = diff(constr_By,2);

odes = [motion_1x; motion_1y; motion_1m; motion_2x; motion_2y; motion_2m; constr_Ax; constr_Ay; constr_Bx; constr_By]

%%b.
cond = [theta1(0) == sym(pi)/2, theta2(0) == sym(pi)/2, dtheta1(0) == 0, dtheta2(0) == 0];
% sol_b(t) = dsolve(odes, cond)

Thanks in advance for trying to help me with this, it is very much appreciated! Cheers

Update: Found out that I don't need to solve the differential equations using dsolve since I want to know the 2nd derivative of a variable instead of the other way around. Thus, I can treat the equations just as normal equations! Silly me

Answer 1

[Torsten]

constr_Ax, constr_Ay, constr_Bx and constr_By give you x1, x2, y1 and y2 depending on theta1 and theta2 (linear system of equations).

Thus you only need to solve motion1m and motion2m together with your initial conditions, and this should make no further problems.

Best wishes

Torsten.

[Vianne Heusdens]

Thanks Torsten! I get your solution, the problem is that filling in the initial conditions by assigning the values to the parameters at time 0 doesn't work as Matlab then sees them as constants and therefore all derivatives become 0. Someone just pointed out to me to substitute all derivatives with some independent parameters (diff(theta) becomes dtheta par example) and then assign values. I hope that, in combination with your solution, will work out!

[Vianne Heusdens]

It did :D