How to use ode45 to solve coupled second order differential equations

Hey, I am conducting a frequency analysis on a system of coupled, second order differential equations. I have come up with this code, but I know that it does not solve my system correctly. With a initial condition of (only) y(0) = 0.2, I know from both physical experimentation and from intuition just from looking at the maths, that both of the systems should oscillate and interact with each other. This is not what is seen here, the y graph clearly shows a purely sinusoidal motion, while the other shows the beating behaviour I am looking for. From this, I know that desolve is not solving my problem correctly.

I have seen posts about solving second order DEs with ode45 but I have no idea what is going on, I don't understand the code at all. So my question is, how could I adapt my code to use ode45 to obtain a solution to these two coupled second order differential equations? Would the fourier() function still work, or would I use fft() when I have obtained numerical data?

clear
syms t omega T
syms y(t) z(t)
syms k m w J I
Dy = diff(y);
Dz = diff(z);
k = 4.5682;
m = 0.309 + 0.032*4; %base mass plus extra experimental mass
w = 0.03612;
J = 0.0327;
I = 5.59*10^-3 + 0.032*4*0.155^2; %base rotational inertia plus extra experimental rotational inertia
odes = [diff(y,2) == (-k/m)*y + (w/m)*z; diff(z,2) == (-J/I)*z + (w/I)*y];
cond1 = y(0) == 0.2;
cond2 = Dy(0) == 0;
cond3 = z(0) == 0;
cond4 = Dz(0) == 0;
conds = [cond1 cond2 cond3 cond4];
qwe = dsolve(odes, conds);
wilberforce_y = simplify(qwe.y, 500)
%Fy = simplify(fourier(wilberforce_y), 500)
wilberforce_z = simplify(qwe.z, 500)
%Fz = simplify(fourier(wilberforce_z), 500)
figure
fplot(wilberforce_y, [0 50])
ylim([-1 1])
grid

figure
fplot(wilberforce_z, [0 50])
ylim([-1 1])
grid
% I had some frequency analysis stuff down below, but it is not necessary
% for the question.

[Steven Lord]

Can you post the mathematical form of the system of ODEs that you're trying to solve? This will help us check the problem from the beginning, making sure that the contents of the odes variable matches the equations you're trying to solve. [If I had a dollar for every time something that I was certain wasn't the problem turned out to be the problem ...]

[Paul]

Amplifying on Steven's comment, dsolve() is correctly solving the differential equations with initial conditions as defined in the problem

clear
syms t omega T
syms y(t) z(t)
syms k m w J I
Dy = diff(y);
Dz = diff(z);
k = 4.5682;
m = 0.309 + 0.032*4; %base mass plus extra experimental mass
w = 0.03612;
J = 0.0327;
I = 5.59*10^-3 + 0.032*4*0.155^2; %base rotational inertia plus extra experimental rotational inertia
odes = [diff(y,2) == (-k/m)*y + (w/m)*z; diff(z,2) == (-J/I)*z + (w/I)*y];
cond1 = y(0) == 0.2;
cond2 = Dy(0) == 0;
cond3 = z(0) == 0;
cond4 = Dz(0) == 0;
conds = [cond1 cond2 cond3 cond4];
qwe = dsolve(odes, conds);
odes = odes(t); % make odes into an array, instead of a function

% Show that the solution satisfies the differential equations
[simplify(diff(qwe.y,t,2) - subs(rhs(odes(1)),[y z],[qwe.y qwe.z])), 
simplify(diff(qwe.z,t,2) - subs(rhs(odes(2)),[y z],[qwe.y qwe.z]))]

ans = 

% Show that the solution satisfes the initial conditions
simplify(subs([qwe.y qwe.z diff(qwe.y) diff(qwe.z)],t,0))

ans = 

[Star Strider]

Reference — Fourier transform of symbolic function

[Bjorn Gustavsson]

The ideal Wilberforce pendulum (i.e. without any damping-factors) leads to two coupled second-order ODEs. They have rather standard harmonic solutions, but rather tedious initial condition set-up: Wilberforce pendulum paper

[Luke Weston]

Yeah, I agree @Bjorn Gustavsson, for context I am trying to improve upon current literature on the Wilberforce pendulum for a u18 award. It would be really helpful if you know anything about the Wilberforce pendulum that would be cool and interesting for me to cover, eg an area of the Wilberforce pendulum that not many papers have covered.

As of now I have done some experimental work, obtained some experimental position/angle data, and after taking the fourier transform of that data, come up with frequency data for each of those tests. I have compared my theoretical and experimental graphs of (mass to inertia to ratio of frequencies) in order to investigate where the point of resonance is for each mass/inertia test. I have also applied my findings to real life by detailing how my findings could be used by engineers to move towards/away from resonance, and what resonance means physically.

Thanks!

[Bjorn Gustavsson]

You could for example take a look at the energy in the two modes (potential and kinetic energy of y and z respectively and look at how the effect of the coupling constant impact the variation of them). You could also take a look at how damping-factors impact the solutions (the original equations are ideal with no dissipation, but where I live everything comes to rest eventually). For that it seems that the symbolic toolbox struggles a bit, dsolve produces solutions with "Root(..." that are difficult to push further, it is possible to calculate the eigenvectors and eigen-values to the 4-by-4 fundamental matrix of the system of equations - but the expressions in them are "rather long". But you might still get some interesting behaviour with numerical solutions.

[Paul]

Getting back to @Steven Lord's comment, are these the correct equations for the Wilberforce pendulum? Should the last terms on the rhs of the equations be:

+ (-w/m)*z*y 

adn

 + (-w/I)*z*y

[Bjorn Gustavsson]

@Paul and @Luke Weston: Check equations 2 and 3 in the paper I linked to above. There seems to be a sign-error in the ODEs above, but the coupling-terms are linear.

[Paul]

Interesting. This reference has the nonlinear terms I cited, but that seems to be an error when they pulled the equations from the the exact reference that you (@Bjorn Gustavsson) cite. I agree that the sign in those terms in the ODEs in the Question should be reversed. At least I had that part correct ;).

[Bjorn Gustavsson]

The way I understood the system is that when a spring (coil) is twisted it also changes its length, and for that type of system a simlpe linear coupling seems to make sense. Perhaps the sign of the coupling-constant depends on the direction of the angular variable - with or against the direction of the spiral.

Answer 1

[Bjorn Gustavsson]

It might be that your expectations on the modification of the beat in y, it might be something else. (When I converted the dsolve solutions to matlab-functions and tried to plot them between 0 and 100 plot warned me about y being complex.) When looking at the real components of the solutions they look like you describe. But when I look at their fourier-transforms it is clear that there is a small coupling:

Wy = matlabFunction(qwe.y);
Wz = matlabFunction(qwe.z);
t = linspace(0,100,5001);
k = 4.5682;
m = 0.309 + 0.032*4; %base mass plus extra experimental mass
w = 0.03612;
J = 0.0327;
I = 5.59*10^-3 + 0.032*4*0.155^2;
subplot(2,2,1)
plot(t,real(Wy(I,J,k,m,t,w)))
subplot(2,2,3)
plot(t,real(Wz(I,J,k,m,t,w)))
subplot(2,2,2)
semilogy(abs(fft(real(Wy(I,J,k,m,t,w)))))
axis([0 100 1 300])
subplot(2,2,4)                           
semilogy(abs(fft(real(Wz(I,J,k,m,t,w)))))
axis([0 100 1 300])

Maybe you can look at the difference between this y and a solution twithout coupling (for that you might need to slightly adjust the spring-constant to get exactly the same frequency?).

HTH

[Luke Weston]

Hey, thanks for the suggestion, however when I run that code, it comes up with the error "Too many input arguments" for the lines containing the bit of code "real(Wy(I,J,k,m,t,w))". Any idea how to fix this?

You can see the code I have integrated with yours below.

clear
syms t omega T
syms y(t) z(t)
syms k m w J I
Dy = diff(y);
Dz = diff(z);
k = 4.5682;
m = 0.309 + 0.032*4;
w = 0.03612;
J = 0.0327;
I = 5.59*10^-3 + 0.032*4*0.155^2;
odes = [diff(y,2) == (-k/m)*y + (w/m)*z; diff(z,2) == (-J/I)*z + (w/I)*y];
cond1 = y(0) == 0.2;
cond2 = Dy(0) == 0;
cond3 = z(0) == 0;
cond4 = Dz(0) == 0;
conds = [cond1 cond2 cond3 cond4];
qwe = dsolve(odes, conds);
Wy = matlabFunction(qwe.y);
Wz = matlabFunction(qwe.z);
t = linspace(0,100,5001);
k = 4.5682;
m = 0.309 + 0.032*4; %base mass plus extra experimental mass
w = 0.03612;
J = 0.0327;
I = 5.59*10^-3 + 0.032*4*0.155^2;
subplot(2,2,1)
plot(t,real(Wy(I,J,k,m,t,w)))
subplot(2,2,3)
plot(t,real(Wz(I,J,k,m,t,w)))
subplot(2,2,2)
semilogy(abs(fft(real(Wy(I,J,k,m,t,w)))))
axis([0 100 1 300])
subplot(2,2,4)                           
semilogy(abs(fft(real(Wz(I,J,k,m,t,w)))))
axis([0 100 1 300])

[Bjorn Gustavsson]

My bad. I messed up with assigning numerical values to k, m, I, J and w before the call to dsolve - then the solution will be for those specific values. If you modify the code to:

clear
syms t omega T
syms y(t) z(t)
syms k m w J I
Dy = diff(y);
Dz = diff(z);
odes = [diff(y,2) == (-k/m)*y + (w/m)*z; diff(z,2) == (-J/I)*z + (w/I)*y];
cond1 = y(0) == 0.2;
cond2 = Dy(0) == 0;
cond3 = z(0) == 0;
cond4 = Dz(0) == 0;
conds = [cond1 cond2 cond3 cond4];
qwe = dsolve(odes, conds);
Wy = matlabFunction(qwe.y);
Wz = matlabFunction(qwe.z);
t = linspace(0,100,5001);
k = 4.5682;
m = 0.309 + 0.032*4; %base mass plus extra experimental mass
w = 0.03612;
J = 0.0327;
I = 5.59*10^-3 + 0.032*4*0.155^2;
subplot(2,2,1)
plot(t,real(Wy(I,J,k,m,t,w)))
subplot(2,2,3)
plot(t,real(Wz(I,J,k,m,t,w)))
subplot(2,2,2)
semilogy(abs(fft(real(Wy(I,J,k,m,t,w)))))
axis([0 100 1 300])
subplot(2,2,4)                           
semilogy(abs(fft(real(Wz(I,J,k,m,t,w)))))
axis([0 100 1 300])

Then it works.

To see that the coupling is "strong" for z but weak for y you can also look at the coefficients for each term:

[k/m w/m J/I w/I]

for z they are rather equal for y quite different.

Another idea is to compare the energies of the different modes - if the energy in the z-mode is way smaller than the energy of y then the perturbation of y will have to be smallish...

HTH