# How to use ode45 to solve coupled second order differential equations

2 views (last 30 days)
Luke Weston on 3 Nov 2021
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.
Bjorn Gustavsson on 5 Nov 2021
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.

Bjorn Gustavsson on 3 Nov 2021
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
##### 2 CommentsShowHide 1 older comment
Bjorn Gustavsson on 3 Nov 2021
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