Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Help with Matlab solving second order differential equations

  1. Oct 19, 2010 #1
    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 '),
    plot(t, x(:, 2), 'b')
    title('Nucleus position vs time'), xlabel('Time '),

    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!
  2. jcsd
  3. Oct 25, 2010 #2

    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.
  4. Oct 25, 2010 #3
    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!
  5. Oct 25, 2010 #4
    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, lets 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);
  6. Oct 25, 2010 #5
    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), becuase 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);
  7. Oct 25, 2010 #6
    What is the value of n? It is likely that n>m
  8. Oct 25, 2010 #7
    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..
  9. Oct 25, 2010 #8
    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).
  10. Oct 25, 2010 #9
    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.
  11. Oct 25, 2010 #10
    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.
  12. Oct 25, 2010 #11
    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!
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook