How to use ode45 to solve coupled second order differential equations

40 visualizaciones (últimos 30 días)
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.
  10 comentarios
Paul
Paul el 4 de Nov. de 2021
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
Bjorn Gustavsson el 5 de Nov. de 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.

Iniciar sesión para comentar.

Respuesta aceptada

Bjorn Gustavsson
Bjorn Gustavsson el 3 de Nov. de 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 comentarios
Luke Weston
Luke Weston el 3 de Nov. de 2021
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
Bjorn Gustavsson el 3 de Nov. de 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

Iniciar sesión para comentar.

Más respuestas (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by