help-octave
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: Trying to solve du/dt = i*u


From: John B. Thoo
Subject: Re: Trying to solve du/dt = i*u
Date: Wed, 13 May 2009 04:54:32 -0700


On May 1, 2009, at 4:09 PM, Carlo de Falco wrote:

On 1 May 2009, at 19:37, John B. Thoo wrote:

That's very neat.  I tried it and it works! :-)  But now how do I
solve a system like this:

udot(1) = -i*u(3)*u(2)' - i*u(2)*u(1)';
udot(2) = -2*i*u(3)*u(1)' - i*u(1)*u(1);
udot(3) = -3*i*u(2)*u(1);

(I'm using  '  for complex conjugate.  Is that correct?)

I don't see how to index the real and imaginary parts of the
functions u(1), u(2), and u(3).

Thank again.

---John.


try the following:

----------8<------------
function ydot = fun(y, t)

u = y(1:3) + i * y(4:6);

ydot(1) = real(-i*u(3)*u(2)' - i*u(2)*u(1)');
ydot(2) = real(-2*i*u(3)*u(1)' - i*u(1)*u(1));
ydot(3) = real(-3*i*u(2)*u(1));

ydot(4) = imag(-i*u(3)*u(2)' - i*u(2)*u(1)');
ydot(5) = imag (-2*i*u(3)*u(1)' - i*u(1)*u(1));
ydot(6) = imag (-3*i*u(2)*u(1));

endfunction

[x, istate, msg] = lsode ('fun', ones(6,1), linspace(0,10,100))
plot(x(:,1), x(:,4), 'r', x(:,2), x(:,5), 'b', (x(:,3), x(:,6), 'k')

----------8<------------
c.

Thanks, Carlo. That worked (no errors). But now I'm having trouble generalizing.

I'm trying to solve  u_t = -(u^2/2)_x  by solving the system of ODEs

                          oo
     d uhat_k        ik   ---
(*)  --------  =  - ----  \    uhat_{k-m} uhat_m
        dt           2    /__
                        m = -oo

[ \frac{d \hat{u}_k}{dt}
   = -\frac{ik}{2} \sum_{m=-\infty}^\infty \hat{u}_{k-m} \hat{u}_m ]

where the  uhat  are functions of  t.  Then

             oo
             ---
  u(x,t)  =  \    uhat_k * e^{ikx}
             /__
            k = -oo

[ u(x,t) = \sum_{k=-\infty}^\infty \hat{u}_k e^{ikx} ].

The system Carlo helped me with (above) is explicitly for k = 1..3, but I don't know how to write the summation in (*) for large k = K (truncating the series by setting uhat = 0 for k > K).

I tried using shift (x, b) but that didn't seem to work. Not being very familiar with Octave, I'm at a loss to know how to proceed. Help?

Thanks very much.

---John.


reply via email to

[Prev in Thread] Current Thread [Next in Thread]