Trying to solve du/dt = i*u

[John B. Thoo]

On May 13, 2009, at 4:54 AM, John B. Thoo wrote:

> 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.

Hi.  I wrote this script to implement Carlo's code in the hope that I  
would be able to generalize it:

%----------begin script----------
% script burgers_s ()

1;   % prevent Octave from thinking this is a function file


% define the space mesh
t0 = 0;   % initial time
tf = 10;   % final time
a = -1;   % left end point
b = 1;   % right end point
M = 100;   % number of time intervals
N = 100;   % number of space intervals
dt = (tf - t0)/M;
t = t0:dt:tf;
dx = (b - a)/N;
x = a:dx:b;

% define the ODE
function ydot = f (y, t)

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

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

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

endfunction

% define initial conditions
yy = (x <= 0);
fftyy = fft (yy);
y0 = horzcat (real (fftyy(2:4)), imag (fftyy(2:4)));
t = linspace (t0, tf, M+1)';

% solve the ODE
[w, istate, msg] = lsode ("f", y0, t);

% recover  uhat
uhat = zeros (3, M+1);
for j = 1:3
   uhat(j, :) = w(:, j)+i.*w(:, j+3);
endfor

% define  u
u = zeros (N+1, M+1);
EXP = zeros (N+1, 3);

for j = 1:3
   EXP(:, j) = exp (i.*j.*x);
endfor

for j = 1:M+1
   u(:, j) = EXP*uhat(:, j);
endfor
%----------end script----------


The script *seems* to work correct.  (If not, please let me know.)   
My problem now is that when I replace the function ydot with this:

%----------begin replacement function----------
function ydot = f (y, t)

K = 3;   % truncation index

uhat = y(1:2*K) + i*y(2*K+1:2*K+2*K);
uhat(K+1:K+K) = 0;

for k=1:K
   for j=1:K
     xxdot(k) = real (-i*k*(2*uhat(k+K+2-j)*uhat(k+2-j)'));
   endfor
   for j=1:k+1
     zzdot(k) = real (-i*k*(uhat(k+2-j)*uhat(j)'));
   endfor
   ydot(k) = xxdot(k) + zzdot(k);
endfor

for k=1:K
   for j=1:K
     xxdot(k+K) = imag (-i*k*(2*uhat(k+K+2-j)*uhat(k+2-j)'));
   endfor
   for j=1:k+1
     zzdot(k+K) = imag (-i*k*(uhat(k+2-j)*uhat(j)'));
   endfor
   ydot(k+K) = xxdot(k) + zzdot(k);
endfor

endfunction
%----------end replacement function----------

I get the following errors:

octave-3.0.5:68> burgers_s
error: invalid vector index = 12
error: evaluating binary operator `*' near line 23, column 20
error: evaluating binary operator `+' near line 23, column 17
error: evaluating assignment expression near line 23, column 6
error: called from `f'
error: lsode: evaluation of user-supplied function failed
error: lsode: inconsistent sizes for state and derivative vectors
error: near line 55 of file `/Users/jbthoo/Desktop/NumMethPractice/ 
burgers_s.m'
octave-3.0.5:68>


What am I doing wrongly?

Thanks.

---John.