numerical settings to improve/stabilize pdepe solution

Hi,

I am using MATLAB (R2017a) pdepe to solve a PDE with third-order in 1-D space and first-order in time. Please see attached for the specific form of the PDE. The coefficients (F1, F2, F3 and C1, C2, C4) are given and in general complex quantities as a function of time (please see attached). Because the coefficients depend on time, I solve the set of equations in a piecewise way, i.e., I use pdepe to solve in the interval [t0,t1], withing which the coefficients are constant (but a complex number). Then I record the ouput [u1,u2,u3] at t1, serving as input to [t1,t2], and proceed with pdepe.

Without adding the higher order spatial derivatives, things work fine, i.e., pdepe can successfully finish at all t. However, when adding the second and third derivatives, at a certain t (not from very beginning), the numerical computing is stuck for a while and then shows the following warning message:

"Warning: Failure at t=4.060671e-03. Unable to meet integration tolerances without

reducing the step size below the smallest value allowed (1.387779e-17) at time t."

The computing efficiency after the warning message becomes much slower and I am not sure if pdepe gives the correct result.

I have three questions specific to the numerical setting for the set of equations:

A) Are the three arrays c, b, and s set correctly or in a favorable way for pdepe?

A-1) I ever tried to replace DuDx(1) and DuDx(2) with u(2) and u(3) and leave other settings unchanged. Then pdepe responds differently: using u(2) and u(3) appears to propagate the system longer than using DuDx(1) and DuDx(2), although both attempts fail to finish the solution. Is there any difference in using the two formats?

A-2) It seems that the second and third equations are much more simple than the first equation and there is no time derivative. Is it suggested to solve the two simpler equations separately? (if so, how?)

B) Are the initial conditions (ICs) and boundary conditions (BCs) set correctly or in a favorable way for pdepe? I ask this because for ICs of u2 and u3, I simply take first and second spatial derivative to u1. Will that lead to undesired numerical effect?

C) Will the choice of xmesh and tspan (especially their ratio) affect the (numerical) stability of spatial-temproal discretization? When I try to vary the number of meshes/grids in x and t, in some situations pdepe can propagate longer than other situations, which was not expected for me. I was wondering if there are some general advices for numerical setting for PDE with higher order of spatial derivatives.

Any help, hints or suggestions are appreciated...

--------------------- Edited on June 19, 2019 ---------------------

Below are the scripts I used to solve the set of equations. The first one works better than the second one, but both fail to finish the calculation. The input data, which provide the time-varying complex coefficients (F1,F2,F3 and C1,C2,C4), are in the attachment. I still wish to seek for any hints or comments on the aforementioned issues.

[scripts removed]

Plus one interesting observation:

D) When I change the order of elements of c,b,s arrays, I find that pdepe gives different results and in some case it even shows an error message:

"Spatial discretization has failed. Discretization supports only parabolic and elliptic equations, with flux term involving spatial derivative."

Still any help, hints or suggestions are appreciated.

----------------------------------------------------------------------------

--------------------- Edited on June 24, 2019 ---------------------

According to Bill's suggestion, I reduced the three equations to two by removing an irrelevant one. Below is the updated code; the updated document includes the equations with initial and boundary conditions. With very short integration distance, e.g., z from 0 to 1, it works. However when increasing the integration range, e.g., z from 0 to 5 (with 1000 grid points), it does not successfully finish and shows the warning message

"Warning: Failure at t=2.738149e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (7.105427e-15) at time t."

clear all
format long

mc2=0.511;                                      % electron rest mass energy, MeV
deg2rad=pi/180;

k_0_tmp=load('root_k_vs_z.txt');
k_0=k_0_tmp(:,2)+1i*k_0_tmp(:,3);               % load from external file
k_0=k_0.';                                      % use .' for non-conjugate transpose

K_0=3.5;                                        % undulator parameter
Theta_R=-90*deg2rad;                            % ponderomotive phase, -90 deg for untapered case
rho=1.57e-3;
A_sat=1;

z_start=0; z_end=15; z_sat=6.0; z_num=500;
z_array=linspace(z_start,z_end,z_num);
[~,ind]=min(abs(z_array-z_sat));

fR=1-2*rho*A_sat*cos(Theta_R)*(z_array-z_sat)-rho*(cos(Theta_R))^2*(z_array-z_sat).^2;  % KMR, constant \Theta_R
fR(1:ind)=ones(1,ind);
fB=sqrt(2)/K_0*sqrt((1+0.5*K_0^2)*fR.^2-1);
etaR=(fR-1)/rho;

tmp01=3*k_0.^2-4*etaR.*k_0+2*etaR.^2;
C1=2*k_0.*fB*cos(Theta_R)./tmp01;
C2=2*1i*(2*k_0.^2+abs(k_0).^2-2*k_0.*etaR)./tmp01;
C4=-2*1i*k_0./tmp01;

% read F_n from external file
Re_F_n_tmp=load('Re_n_vs_z.txt');
Im_F_n_tmp=load('Im_n_vs_z.txt');
F1_array=Re_F_n_tmp(:,3)+1i*Im_F_n_tmp(:,3);
F2_array=Re_F_n_tmp(:,4)+1i*Im_F_n_tmp(:,4);
F3_array=Re_F_n_tmp(:,5)+1i*Im_F_n_tmp(:,5);

% ============================= MAIN ==================================== %
m_order=0; 
tau_num=1000; tau_scale=15; T0=1;
tau_array=linspace(-tau_scale*T0,tau_scale*T0,tau_num); 

z_start_CGLE=0; z_end_CGLE=5; z_num_per_seg=1000;
Z_array=linspace(z_start_CGLE,z_end_CGLE,z_num_per_seg);

sol=pdepe(m_order,@eqn_CGLE_rd_stf,@initial_CGLE_rd_stf,@bc_CGLE_rd_stf,tau_array,Z_array,[],z_array,k_0,F1_array,F2_array,F3_array,C1,C2,C4);
    
U_fun=sol(:,:,1);
ini_fun(1,:)=sol(end,:,1);
ini_fun(2,:)=sol(end,:,2);

Int_fun=abs(U_fun).^2;
figure(4); surf(tau_array,Z_array,Int_fun); xlabel('\tau'); ylabel('z'); 
title('Intensity (arb. units)'); shading interp; colorbar; view(0,90);
% ======================================================================= %

function [c,b,s]=eqn_CGLE_rd_stf(tau,z,u,DuDx,z_array,k_0,F1_array,F2_array,F3_array,C1,C2,C4)
tmp01=interp1(z_array,k_0,z);
tmp02=interp1(z_array,F1_array,z);
tmp03=interp1(z_array,F2_array,z);
tmp04=interp1(z_array,F3_array,z);
tmp05=interp1(z_array,C1,z);
tmp06=interp1(z_array,C2,z);
tmp07=interp1(z_array,C4,z);

k0_R=real(tmp01); 
k0_I=imag(tmp01);

c=[1;0];
b=[(k0_R*tmp02*u(1)+k0_R*tmp03*DuDx(1)+k0_R*tmp04*DuDx(2));-u(1)];
s=[1i*tmp01*u(1)+tmp05*abs(u(1))*u(1)+tmp06*abs(u(1))^2*u(1)+tmp07*abs(u(1))^4*u(1);u(2)];
end

function [pl,ql,pr,qr]=bc_CGLE_rd_stf(taul,ul,taur,ur,z,z_array,k_0,F1_array,F2_array,F3_array,C1,C2,C4)
pl=[ul(1);ul(2)];
ql=[0;0];
pr=[ur(1);ur(2)];
qr=[0;0];
end

function value=initial_CGLE_rd_stf(tau,z_array,k_0,F1_array,F2_array,F3_array,C1,C2,C4)
sig_tau=1;
tmp00=exp(-tau^2/(2*sig_tau^2));
tmp01=tmp00/sqrt(2*pi)/sig_tau;
tmp02=tmp00*tau/sqrt(2*pi)/sig_tau^3;
value=0.1*[tmp01;tmp02];
end

Any hints or further suggestions are appreciated...

----------------------------------------------------------------------------

[Bill Greene]

It is unclear to me why you are stopping and restarting pdepe? pdepe can handle terms that are functions of time in a straightforward way.

[Cheng-Ying Tsai]

Hi Bill,

Thanks for the comment. When solving in a straightfforward way, where the coefficients are evaluated by interpolation, it stops at very early stage. Then I ever tried to fix the coefficients to constant values and still solve in a straightfforward way, it works. I thought that pdepe might jump back and forth in a small range of time, where it will sample different values of the coefficients, and could not easily find a stable/converged solution. Therefore I rewrite the script to solve the set of equations in a piecewise way, and it works fine but only for u(2) and u(3) [or DuDx(1) and DuDx(2)] not included.

Any hints, comments or suggestions for other parts?

Thanks!

[Bill Greene]

Your code doesn't agree with your mathematical description and this makes it difficult to understand and, I believe, has also lead to errors. One thing I noticed immediately is that the code models the domain as zero to some length and the mathematical description shows the domain as -L to +L.

If you plot your initial conditions, as I did, you will see that they are obviously incorrect.

In my experience, the only way to debug problems of this complexity is to systematically check each small piece of the problem individually.

[Cheng-Ying Tsai]

Hi Bill,

Thanks for your comments. The domain I am modeling refers to "tau_array", which is from -15 to +15. The temporal evolution is from "z_start_CGLE" to "z_end_CGLE" in the code, which is from 0 to 15. The initial condition is assumed Gaussian. I plot it using the command "plot(tau_array,ini_fun(1,:))" and it works. Could you elaborate where the erroes are?

Thanks very much for your interest and help!

Cheng-Ying

[Bill Greene]

Oh, sorry, I was confused by the inconsistencies in the notation between your document and your code.

[Cheng-Ying Tsai]

Hi Bill, sorry for confusion in the notation used in the script and mathematical formulas. I appreciate if any your (or other experts') comments or hints for the above questions (A) to (D). I've tried but have no idea why pdepe does not work throughout.

[Bill Greene]

Regarding your document, I don't see a reason for the third equation so I suggest you remove that. Other than that change, it looks OK to me.

I see no reason for breaking the time interval into segments. So I suggest you abandon that approach.

At this point, I suggest you make the changes I suggest to the doc and then implement a matlab version that corresponds as closely as possible to that document both in content and notation. Edit your post to include that version of the doc and code.

[Bill Greene]

I don't understand your initial conditions. tmp02 doesn't look like the derivative of tmp01 wrt to tau to me.

There appears to be a small error in z_array. You calculate it as:

z_start=0; z_end=15; z_sat=6.0; z_num=500;
z_array=linspace(z_start,z_end,z_num);

But from your data files (e.g. root_k_vs_z.txt) it looks like the data actually starts at z=.03. Shouldn't it be:

z_array=k_0_tmp(:,1)';

This change will require interpolation as, e.g.

tmp01=interp1(z_array,k_0,t,'linear','extrap');

At this point, I'm simply trying to eliminate all the errors I can find in the code to see if they affect the numerical instability that develops. These types of problems can be very sensitive to initial and boundary conditions.

[Cheng-Ying Tsai]

Hi Bill,

Thanks for the comment; I will test it tonight according to your advice, especially the interpolation part. Let me quickly respond to an issue you raised:

about initial condition: do you mean a missing negative sign? Indeed I missed it (I will append it).

Thanks.

[Bill Greene]

After looking at the initial conditions a second time, the only issue I see is the negative sign. And, in fact, I don't think it matters much to the solution because I believe, internally, the code is modifying the initial conditions to obtain values that satisfy the PDE at the z=0.

[Cheng-Ying Tsai]

Hi Bill,

Thank you very much for your comments and pointing out some bug in my code. After adding a negative sign to initial condition u_2(t=0,\tau) and linear extrapolating the coefficients, it still gets stuck at some intermediate step. Eventually the code displays a warning message and fails to finish the calculation. It could get stuck around z=1.2, after I make the above two modifications to the "piecewise" version of the code.

I think I may consult more specific numerical algorithms to this particular equation, to see if there is some way to work around this. But, I still seek for (if any) possible hints/comments/suggestions from you or other experts. Thanks!

[Bill Greene]

Unfortunately I have not found a way to stabilize this equation.

The reason ode15s fails is the spurious oscillations typical of what

one would see when trying to solve the convection-diffusion equation

with pdepe. The workaround in the convection-diffusion problem is to

add some additional, artificial diffusion. For the current equation, the

and terms are actually negative over parts of their

time history. This tends to destabilize the solution. I don't believe pdepe

is the right solver for this equation. And am unaware of any PDE solvers

that can handle complex coefficients and also cope with the hyperbolic

character of this equation.

[Cheng-Ying Tsai]

Hi Bill,

Thank you for the comment. I did what you mentioned in adding a small artiifcial diffusion term, as you suggested in other/earlier similar posts, but it fails. I agree with you that and (ps. I just know it supports LaTeX :) play a key role to the dynamical behavior of the solution (e.g., as you pointed out, the negative part). I am searching for other possible numerical approach to solve this type of equation. Again, thank you very much for your interest, time, and help in this situation.

Cheng-Ying