Help with Matlab solving second order differential equations

[Davide86]

I am a Matlab rookie. I need to solve numerically the following second order differential equations

d^2x/dt^2 + w0_(el) * x = e/m_e * E - K3/m_e * x *y;
d^2y/dt^2 + w0_(v) * y = - K_3/2M * x^2;

I have started to deal with only the harmonic part of the problem. So I tried to solve

d^2x/dt^2 + w0_(el) * x
d^2y/dt^2 + w0_(v) * y

with the following program

t = 0:10^(-15):0.5*10^(-9);
x0 = zeros(1,2);
x0(1) = input('Insert the initial value of x');
x0(2) = input('Insert the initial value of dx/dt');
[t, x] = ode45(@harmonic, t, x0);
plot(t, x(:, 1), 'g')
title('Electronic position vs time'), xlabel('Time '),
ylabel('Position')
hold
figure
plot(t, x(:, 2), 'b')
title('Nucleus position vs time'), xlabel('Time '),
ylabel('Position')

where the function "harmonic" is

function ydot = harmonic(t, x)
ydot = zeros(2,1);
w_e = 10^30;
w_v = 10^24;
ydot(1) = x(3);
ydot(2) = x(4);
ydot(3) = -w_e*x(1);
ydot(4) = -w_v*x(2);

The outcome is a damping oscillation behaviour for x, which is of course meaningless because there aren't damping terms in the harmonic equations. So I am pretty desperate, if someone can help me I wuold be very glad.
Thank you!

 

[MasterX]

Hi,

I believe that you have mixed your variables in the ydot function. Try the following ydot function

function ydot = harmonic(t, x)
ydot = zeros(2,1);
w_e = 10^30;
ydot(1) = x(2);
ydot(2) = -w_e*x(1);

Write a similar function for the w_v.

 

[Davide86]

Thank you for your help MasterX! What y suggested works, but I need a single expression for both the oscillators. In this way I can add the coupling terms. I found out that the main problem was the solver: the ode45 is not suitable for stiff problems, like the one I am considering.
Thanks again!

 

[MasterX]

Actually, ode45 would do a pretty good job with stiff systems, as it is a higher order method.
As you have written the ydot function would not work. You have mixed the variables

So, let's assume that x(1) is x, x(2) is dx/dt, x(3) is y and x(4) is dy/dt. Note, that you must a similar assignment for the ydot variable.

This is what the ydot function should be:

function ydot = harmonic(t, x)
ydot = zeros(4,1);
w_e = 10^30;
w_v = 10^24;
ydot(1) = x(2);
ydot(2) = -w_e*x(1);
ydot(3) = x(4);
ydot(4) = -w_v*x(3);

 

[Davide86]

I have tried your code: it works. Now I have added the following part, which should provide the Fourier transform od the oscillations. The problem is that I can see two sharp peaks at 10^14 Hz for x(1) oscillations and 10^12 Hz for nuclei: they are shifted by one rder of magnitude in comparison with the input data. I mean, they should and must be at 10^15 and 10^13 Hz (I changed w_v from 10^24 to 10^26), because in the harmonic oscillator equation the frequency are squared and the input values are 10^30 and 10^15. Thank you again!

e = x(:, 1); % electrons' positions
p = x(:, 3); % phonons' positions
m = length(e);
n = pow2(nextpow2(m));
y = fft(e, n);
q = fft(p, n);
fs = (5*10^(-11)/n)^(-1); % step in the freq domain
f = (0:(n-1))*(fs/n); % Frequency range
y0 = fftshift(y);
q0 = fftshift(q);
f0 = (-n/2:n/2 - 1)*(fs/n);

 

[MasterX]

What is the value of n? It is likely that n>m

 

[Davide86]

it can be n>m. In that case MATLAB adds zero values to the signal's vector, in order to make the length of the vectors match. At least that's what I have understood from the Mathlab help..

 

[MasterX]

That could be a problem. You should decrease the value of nextpow2(m) by one (or more) so that n<=m. Otherwise, you are taken the FFT of a different function than the one you have (unless, you know that as n->oo, x(:,1)->0).
If n<m then you should not use the first n data of x(:,1) array, but instead, you should try to use n data from the entire x(:,1) array. Namely, if you take the first n data, you are ignoring the effect that the last m-n data could have! Thus, try to take n data, but use both ends of the vector x(:,1).As a first attempt, I suggest to set n=m, and then use the nextpow2(m) function (make sure that n<=m).

 

[Davide86]

I have just try to set m=n and the problem has not been solved. Maybe the problem could be in the definition of the frequency axis f0 or, I think most likely, in the edge effects of the FFT: the two peaks that I see are at the edges of the frequency axis.

 

[MasterX]

Well, I need to study a little more about the FFT :) What are fs,f equal to? I have not used them before!
Try something that is extremely easy. Change the time step, so that m=2^(any number you wish). In this way, pow2(nextpow2(m)) is equal to m.

 

[Davide86]

fs is the step in the frequency axis: 5*10^{-11} is the total temporal window (50 ps), n is the number of the time axis. SO fs is the inverse of the time step.
f is the frequency axis, but for the plot I use f0 because it is zero-centered and it shows the two parts (positive and negative) of the FFT.
Thank you again so much for your help, I owe you at least one beer!