Open boundary in wave equation with “applyBoundaryCondition”

I am implementing a simulation of the wave equation using "solvepde” function in MATLAB. the wave is travelling to right from edge E1 and propagate off forever the edge E3 without being disturbed. how can I specify open boundary (full absorbing/nonreflecting boundary) condition with “applyBoundaryCondition” command?

I am aware there are many academic papers discussing nonreflecting/absorbing boundary conditions but most seem to focus on analytic solutions. On open boundary we have:

I cannot figure out how to implement nonreflecting boundaries numerically in my simulation. This is the code of incident wave on edge E1:

fun1=@(location,state) 2*sin(2*pi/T*state.time)

applyBoundaryCondition(model,'dirichlet','Edge',1,'u',fun1)

Does anyone know how can I adjust the “applyBoundaryCondition” to have open boundaries on edge E3?

[Torsten]

As far as I know, you only have the choice between Dirichlet, Neumann and Mixed for setting up boundary conditions. None of these reflects an "Open Boundary".

[Yaser]

as showen as in figure on open boundary(edge E3) we have du/dx=du/dt (if c=1). how does Dirichlet, Neumann conditin reflect this condition in boundary?

[Torsten]

As I said, MATLAB does not offer an appropriate boundary condition type to set du/dt = c*du/dx at a boundary.

Use a tool especially designed to cope with hyperbolic PDEs, e.g. CLAWPACK available under

https://www.clawpack.org/

[Yaser]

According to Generalized Neumann condition

n·(c×∇u) + qu = g → du/dx+qu=g

if q=0 and g= du/dt

then Neumann condition satisfy open boundary condition (c=1)

du/dx= du/dt

and so how we can assign g as du/dt with applyBoundaryCondition

[Torsten]

As far as I know, g can only depend on x, u and t, not du/dt.

[Yaser]

there are sum commans like diff and gradient which in my opinion can reflect the du/dt, but i get some errors as follow:

fun=@(location,state) diff(state.u,state.time)
applyBoundaryCondition(model,'neumann','Edge',3,'g',fun,'q',0);

error: Difference order N must be a positive integer scalar.

Or

fun=@(location,state) gradient(state.u)/gradient(state.time)
applyBoundaryCondition(model,'neumann','Edge',3,'g',fun,'q',0);

error: Solution failed to reach the requested end time.

[Torsten]

If it says in the documentation that g can depend on x,t and u, then you will have to accept this.

[Yaser]

how this error can be correct?

Warning: Failure at t=0.000000e+00. Unable to meet integration tolerances without reducing

the step size below the smallest value allowed (7.905050e-323) at time t.

> In ode15s (line 661)

In pde.EquationModel/solveTimeDependent (line 101)

In pde.PDEModel/solvepde (line 57)

In for_dr (line 44)

Error using pde.EquationModel/solveTimeDependent

Solution failed to reach the requested end time.

Error in pde.PDEModel/solvepde (line 57)

[u,dudt] = self.solveTimeDependent(coefstruct, u0, ut0, tlist, ...

Error in for_dr (line 44)

result = solvepde(model,tlist);

[Torsten]

If the code works if you change the boundary condition at the exit to a supported type, you know what the reason for the time integration failure is ...

[Yaser]

how can i correct this error?

Warning: Failure at t=0.000000e+00. Unable to meet integration tolerances without reducing

the step size below the smallest value allowed (7.905050e-323) at time t.

Solution failed to reach the requested end time.

[Torsten]

I think nobody here has enough experience with the PDE toolbox to tell the reason for the error - even if we knew your code (which is of course a requirement).

[Yaser]

clc
clear all
close all
m=1

m = 1

d=0

d = 0

c=1

c = 1

a=0

a = 0

f=0

f = 0

% ...........
Domain=[3 4 0 0 10 10 0 12 12 0]';
ns=[];
sf=[];
geom  = decsg(Domain)%decsg(Domain,sf,ns)

geom = 7×4

2 2 2 2 0 0 10 10 0 10 10 0 0 12 12 0 12 12 0 0 0 0 0 0 1 1 1 1

model = createpde; % Set model geometry
geometryFromEdges(model,geom)

ans = 
  AnalyticGeometry with properties:

       NumCells: 0
       NumFaces: 1
       NumEdges: 4
    NumVertices: 4
       Vertices: [4×2 double]

generateMesh(model,'Hmax',.5,'Hmin',.2);hold on
pdemesh(model)
pdegplot(model,'EdgeLabels','on'),
specifyCoefficients(model,'m',m,'d',d,'c',c,'a',a,'f',f);
T=2

T = 2

fun1=@(location,state) 1*sin(                 +2*pi/T*state.time)

fun1 = function_handle with value:
    @(location,state)1*sin(+2*pi/T*state.time)

bc=applyBoundaryCondition(model,'dirichlet','Edge',1,'u',fun1,'Vectorized','off');  %
% % % ==============================================================================================
fun4=@(location,state) -gradient(state.u)./gradient(location.x)

fun4 = function_handle with value:
    @(location,state)-gradient(state.u)./gradient(location.x)

applyBoundaryCondition(model,'neumann','Edge',3,'g',fun4,'q',0);
xlabel x
ylabel y

% initial conditions:
u0=0.1

u0 = 0.1000

ut0=0

ut0 = 0

setInitialConditions(model,u0,ut0);
% initial:
tlist=linspace(0,5);
model.SolverOptions.ReportStatistics ='on';
result = solvepde(model,tlist);

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

0 successful steps
615 failed attempts
2462 function evaluations
1 partial derivatives
616 LU decompositions
2461 solutions of linear systems

Error using pde.EquationModel/solveTimeDependent
Solution failed to reach the requested end time.

Error in pde.PDEModel/solvepde (line 57)
[u,dudt] = self.solveTimeDependent(coefstruct, u0, ut0, tlist, ...

u = result.NodalSolution;

[Yaser]

i have sent the cod.

please run it on your system

[Yaser]

Can anybody solve above error?

[Torsten]

Don't set the boundary condition at edge 3 by fun4. Use a usual Neumann condition.

Further, you didn't set boundary conditions at edge2 and edge4. I don't think the PDE toolbox likes this.

[Yaser]

according to matlab help:

Note that applyBoundaryCondition uses the default Neumann boundary condition with g = 0 and q = 0 for equations for which you do not explicitly specify a boundary condition.

So no boundary condition at edge2 and edge4 wont have any problem.

as mentioned earlier we want to apply open/non reflective boundary condition at edge 4.

[Torsten]

as mentioned earlier we want to apply open/non reflective boundary condition at edge 4.

Edge 4 ? I thought it was edge 3.

And as I mentionned earlier: This is not possible with the PDE toolbox.

Or why do you think the time integrator cannot solve your problem ? It's the boundary condition at edge 3 that causes the problems.

[Yaser]

Sory i had mistake. edge 3 is true. we want to apply open/non reflective boundary condition at edge 3.

Further, we dont use Pde toolbox. we are using solvepde that can solve time dependent problem which PDE toolbox cannot. in matlab help we have:

[Torsten]

"solvepde" is part of the PDE toolbox.

Look at the left-hand side of the page

https://de.mathworks.com/help/pde/ug/pde.pdemodel.solvepde.html

to see this.

I do not doubt that the PDE toolbox can solve the wave equation, but not with a non-reflecting boundary condition.

Try a different boundary condition at edge 3 - if "solvepde" works with the changed setting, you know at least the reason why the integrator fails.