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