[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: Trying to solve du/dt = i*u
From: |
Thomas Shores |
Subject: |
Re: Trying to solve du/dt = i*u |
Date: |
Wed, 13 May 2009 22:53:29 -0500 |
User-agent: |
Thunderbird 2.0.0.21 (X11/20090320) |
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')