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

2 views (last 30 days)

Show older comments

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 Comments

Bjorn Gustavsson
on 5 Nov 2021

### Accepted Answer

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 Comments

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

### More Answers (0)

### See Also

### Community Treasure Hunt

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

Start Hunting!