Weird behavior of ODE45?

% Solve with Fourier-Galerkin with odd expansion
%   u_t =  u_x to show the 

clear; close all;

a=20;       %  domain [-a, a]

N = 128;   % Even N: will use odd-expansion 
T = 2;
dt = 0.1;   % how often we check the numerical solutions.
NStep = round(T/dt);

% Initial condition
uIC = @(x) sin(pi*x/10);


dx = 2*a/(N+1); x = -a+ (0:N)*dx; x = x';
u0 = uIC(x);

b = pi/a;


% IC to frequency space
uF = fft(u0)/(N+1);
uF = [uF(N/2+2:end); uF(1:N/2+1)];  % reording for ODE




% Options for ode45
odeTOL = 1e-10;
odeopt=odeset('AbsTol',odeTOL,'RelTol',odeTOL);
ODEPlist = [b, N];


% 1: Solve with ode45 without interuption, plot the solutions
tspan = dt*(0:NStep);
[locT, uF_new] = ode45(@(t,y)semiEqu_linear(t,y,ODEPlist), tspan, uF, odeopt);
figure(1);
for j=1:NStep
    uNumer  = ifft([uF_new(j+1, N/2+1:end), uF_new(j+1,1:N/2)]')*(N+1);
    plot(x, uNumer);
    title(['Solution from ode45 w/o interupton. T=',num2str(locT(j+1))]);
    pause(0.1);
end


% 2: Solve with ode45, but now stop and check (plot) the solution for every
%  time elapse 0.1
for k=1:NStep
    tspan = [dt*(k-1),  dt*k]    
    [locT, uF_new] = ode45(@(t,y)semiEqu_linear(t,y,ODEPlist), tspan, uF, odeopt);
    
    uF = uF_new(end,:); uF = uF';
   
    uNumer  = ifft([uF_new(end, N/2+1:end), uF_new(end,1:N/2)]')*(N+1);
    figure(2);
    plot(x, uNumer);
    title(['T=',num2str(locT(end))]);
    pause(0.1)    
end % for k=1:NStep


% semi-discrete system 
function dy = semiEqu_linear(t,y, ODEPlist)
b = ODEPlist(1);
N = ODEPlist(2);

dy = i*b*y.*((-N/2:N/2)');      % for u_t = u_x
end

Use Fourier-Galerkin to solve u_t = u_x. The function "semiEqu_linear" defines the semi-discrete system, resulted from Galerkin projection.

1) Solve with ode45 at dt*(0:NStep). So no stopping in the middle to check. You can see the solution is a traveling wave (to the right).

2) Still solve with ode45, but now stop at every dt*k to plot it. Basically, everytime, we use ode45 to solve with t0 = (k-1)*dt, and tF = k*dt. This way, the 'final' solution from the last step should be used as the initial condition of the next step. That is why I update 'uF', which is the initial condition in the ode45 call.

Somehow, when we do this way, we don't see a traveling wave. Is there a bug in the above code? Or ode45 and many other solves are designed this way?? I have to agree it is a bit odd to use ode45 this way.

[Torsten]

[locT, uF_new] = ode45(@(t,y)semiEqu(t,y,ODEPlist), tspan, uF, odeopt);

instead of

[locT, uF_new] = ode45(@semiEqu, tspan, uF, odeopt, ODEPlist);

And you are aware that you never save the solution at the output times because you overwrite uF each time ode45 returns ? So if you want to plot uF after the loop, you'll only see the value at t=Nstep*dt.

[Jan]

Please explain, what this means: "the initial condition 'uF' wasn't updated correctly".

We cannot run your code due to the missing data and functions. If you do not explain, what you observe, we do not have a chance to guess the details.

Answer 1

[Jan]

It is still not clear, what this means: "the initial condition 'uF' wasn't updated correctly"

When I run your code, I do see the changing values of uF. I've removed the pause command and inserted:

disp(max(uF))

"we should see a sinusoidal wave propogate to the right."

We do not see this, but a flipping of the sign.

[EDITED after the question was modified]

:-) Now it gets clear. What is the difference between both codes? It is hidden in this line:

uF = uF_new(end,:); uF = uF';

Remember that ' is the abbreviation for ctranspose - the complex conjugate transpose. This modifies the imaginary parts of the initial values in each step. You want a transposition without conjugating instead:

uF = uF_new(end,:); uF = uF.';
%                          ^   .' : transpose, ' : ctranspose !

Or directly:

uF = uF_new(end,:).';

And then both integrations have the same output.

By the way, the most likely unwanted conjugation happens here again:

uNumer  = ifft([uF_new(end, N/2+1:end), uF_new(end,1:N/2)]')*(N+1);
%                                                         ^

With .' the wave is moving to the left.

[Qianyong Chen]

Thanks for answering. I just edited my original post to make it less confusing.

[Qianyong Chen]

You are right! I had the same issue many years ago, and forgot ' means conjugate transpose.

[Jan]

It was a strange idea to let ' be the conjugate transpose when Matlab was created. And now it must be kept for backward compatibility.

See the last paragraph of my answer: Is ' or .' wanted there?

[Qianyong Chen]

Should be. We don't want to change the components to complex conjugate before applying the inverse discrete Fourier transform.