numerical settings to improve/stabilize pdepe solution

17 次查看(过去 30 天)
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...
----------------------------------------------------------------------------
  13 个评论
Bill Greene
Bill Greene 2019-6-29
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
Cheng-Ying Tsai 2019-6-30
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

请先登录,再进行评论。

回答(0 个)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by