Trying to solve du/dt = i*u

John B. Thoo jthoo at yccd.edu
Wed May 13 06:54:32 CDT 2009


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.


More information about the Help-octave mailing list