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

[Thomas Shores]

John B. Thoo wrote:
> On May 13, 2009, at 4:54 AM, John B. Thoo wrote:
>   
> 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.
> _______________________________________________
> Help-octave mailing list
> address@hidden
> https://www-old.cae.wisc.edu/mailman/listinfo/help-octave
>   
I'm not going to dissect your code, but I will include a template script 
that reduces your problem to making a well-formed definition of the 
complex rhs f(z,t).  It's at the bottom of this email, and if you 
remark/deremark the appropriate lines, you'll get each of the first two 
examples you proposed.

However, I have a few comments about your method, which is a separation 
of variables type argument with complex exponentials whose 
time-dependent coefficients are just the Fourier coefficients of the 
solution:

1. There is an immediate problem with truncation m=-M..M of the double 
series in the indices k, m which comes originally from a double series 
in k and ell, with k+ell=m. The problem is that for all indices k except 
k=0, you will be missing some terms with index k-m, which compounds the 
mathematical truncation error.  It also means you have to be careful 
about limits when you calculate the sums involved in the rhs.

2.  There is a question as to why you are using a Fourier method at all 
for Burgers equation.  This nonlinear hyperbolic problem will generate 
shocks in the solution at a minimum time T = -min u_0(x), where u_0(x) = 
u(x,0) is the initial condition and the derivative  u_0(x)'  is 
somewhere negative.  At this time, the problem becomes ill-posed and 
Fourier expansions are not helpful.  It is this very difficulty that is 
the starting point for the modern theory of "numerical methods for 
conservation laws" (see, e.g., LeVeque's monograph of the same name.)

Hope this helps,

Tom Shores

%-------------begin script------------------
% Script for solving dz/dt = f(z,t), z(t0) = z0
% Here t0 is first coordinate of time node vector T
% at which solution is evaluated.
% User defines f as cplxfcn below as well as
% initial condition and time node vector.
% Vectors z, z0, f(z,t) must be row vectors.

%z0 = 1;
z0 = [1+i,1-i,-2]; % initial condition as row vector
%T = linspace(0,2*pi,20);
T = linspace(0,2,100); % time node vector

function w = cplxfcn(z,t)
 %w = i*z;
 w = i*[z(3)*z(2)'-z(2)*z(1)',-2*z(3)*z(1)'-z(1)*z(1),-3*z(2)*z(1)];
end

global n;
n = length(z0);

function y = realfcn(x,t)
global n;

 z = x(1:n) + i*x(n+1:2*n);
 w = cplxfcn(z,t);
 y = [real(w),imag(w)];
end

% Each row of Z is a time slice of solution with time values from T.
[X, istate, msg] = lsode ('realfcn',[real(z0),imag(z0)], T);
Z = X(:,1:n) + i*X(:,n+1:2*n);

% Write your own output/plot routines for the solution matrix Z, e.g.,
%plot(Z(:,1),'r')
plot(Z(:,1), 'r', Z(:,2), 'g', Z(:,3), 'b')