# How to specify non-reflective boundary conditions when NDSolving two coupled first order PDEs?

Due to a lack of fruitful answers after 3 days, this question has been reformulated in a more direct way compared to its first version.

Consider the telegraph equation problem expressed as a system of two coupled PDEs in voltage v and current i:

eq1 = D[v[t, x], x] + D[i[t, x], t] == 0;
eq2 = D[i[t, x], x] + D[v[t, x], t] == 0;
sys = {eq1, eq2};


The initial conditions are a bell shape centered at x=0.5 for the voltage, and zero current i:

shape = D[0.125 Erf[(x - 0.5)/0.125], x];
ics ={v[0, x] == shape, i[0, x]==0};


We now want to find numerically v and i with NDSolve and the Finite Element method, between xMin and xMax, for t between 0 and 2, and for perfectly reflective or perfectly non-reflective boundary conditions.

xMin = 0; xMax = 1;
region = Line[{{xMin}, {xMax}}];


Case 1: perfect reflection on open ends => i=0

bcs1 = {i[t, xMin] == 0, i[t, xMax] == 0};
sol1 = NDSolve[{sys, ics, bcs1} // Flatten, {v, i}, {t,0,2}, {x} \[Element] region][[1]];


which we plot (current in cyan) using

myInterpPlot[interpolatingFunc_, t_, opts___] :=Block[
{grid = ((#[[2]] &) /@ PropertyValue[interpolatingFunc, Grid][[1]]) //Sort,
points},
points = {#, interpolatingFunc[t, #]} & /@ grid;
ListPlot[{points, points}, Joined -> {True, False}, opts,
PlotStyle -> {{Thin, Red}, {Medium, Blue}}, Frame -> True]]

myPlot [sol_,t_]:=Show[
myInterpPlot[sol[[1, 2]], t, PlotRange -> {-1.5, 1.5}],
myInterpPlot[sol[[2, 2]], t, PlotRange -> {-1.5, 1.5},
PlotStyle -> {{Thin, Cyan}, {Medium, Cyan}}],
Plot[shape, {x, xMin, xMax}, PlotStyle -> {Dashed, Blue}]]

Manipulate[myPlot[sol1,t],{t,0,2}]


Case 2: perfect reflection on shorted ends => v=0

bcs2 = {v[t, xMin] == 0, v[t, xMax] == 0};
sol2 = NDSolve[{sys, ics, bcs2} // Flatten, {v, i}, {t,0,2}, {x} \[Element] region][[1]];
Manipulate[myPlot[sol2,t],{t,0,2}]


Case 3 that I don't know how to encode: perfect non-reflective conditions. Physically, this case corresponds to the two equations {eq1,eq2} applying also on the boudaries, since any wave should propagate as if there were no boundaries. How to encode this case keeping the problem in v and i at first order? Can it be done with bcs3=... ? Should one use Neumann values in eq1 and eq2, but what would they mean for first order equations in space?

Case 4: No specified boundary conditions. What does NDSolve do in this case? What are the default boundary conditions used?

sol = NDSolve[{sys, ics} // Flatten, {v, i}, {t,0,2}, {x} \[Element] region][[1]];
Manipulate[myPlot[sol,t],{t,0,2}]


Thank you for your help. Any reference to the doc is welcome. Note that I know how to encode the non reflective case with a second order eq in v and Neumann values. So this is not the answer to the question.

• Do you want to find a solution for non-reflecting boundary conditions? Mar 14, 2019 at 15:40
• Yes Alex, I want to redo with {eq1,eq2} what I have done on eq with the non reflecting boundary conditions (also called perfectly absorbing or impedence matched). I edit the title of my question and my introduction to be clearer. Mar 15, 2019 at 9:35

Two things:

1. your conversion to a system of first order equations does not work quite right.
2. You'd need to take the time derivative of the help variable out - because you are now dealing with first order time derivatives.

Try this:

solb = NDSolveValue[{D[v[t, x], t] == u[t, x],
D[u[t, x], t] ==
D[v[t, x], {x, 2}] + NeumannValue[-u[t, x], x == 0 || x == 1],
v[0, x] == shape, u[0, x] == 0},
v, {t, 0, 1}, {x} \[Element] region]


Update:

What follows is a telegraph equation like I know it. You can still add coefficients (for dt^2, dt and other components). This equation has a DirichletCondition on the left and an absorbing BC on the right; you'd still need to transform that into a system of first order odes and modify the NeumannValue like in the example above.

shape = D[0.125 Erf[(x - 0.5)/0.125], x];
ufun = NDSolveValue[{D[u[t, x], {t, 2}] + D[u[t, x], t] ==
D[u[t, x], {x, 2}] + u[t, x] +
NeumannValue[- Derivative[1, 0][u][t, x], x == 1],
DirichletCondition[u[t, x] == 0, x == 0],
u[0, x] == shape, Derivative[1, 0][u][0, x] == 0},
u, {t, 0, 2}, {x, 0, 1}];


Look at a particular time:

Plot[ufun[0.58, x], {x, 0, 1}, PlotRange -> {-1.2, 1.2}]


Manipulate[
Plot[ufun[t, x], {x, 0, 1}, PlotRange -> {-1.2, 1.2}], {t, 0, 2}]


Concerning your last question: If no boundary conditions are given on a particular part of the boundary then a NeumannValue[0, part] is implicitly used.

• Thank you user 21 Mar 15, 2019 at 23:03
• Thank you user 21, but I don't understand your answer. The question is about encoding correctly in Mathematica the integration of a system of two coupled first order PDEs with absorbing boundary conditions. You propose to use the second order differential equation, which indeed works as I have already pointed out... You have just hidden the time derivative of i in a help variable u. Mar 15, 2019 at 23:13
• Thank you again. What you solve is a PDE which is second order in space and second order in time.Then, because it becomes second order in space, it has a gradient of a gradient and a Neumann value for the gradient can be given. My question about Mathematica (and not about Physics) is how do you tell Mathematica to NDSolve the system {eq1,eq2} of two coupled first order PDEs, with absorbing boudaries, without translating the problem into a second order equation? Is there a Mathematica syntax for encoding the non-reflective conditions when applying NDSolve to {eq1,eq2}? Mar 18, 2019 at 15:26
• And actually, I want non-reflective conditions on both ends (no Dirichlet condition) but it is a detail. Mar 18, 2019 at 15:28
• And concerning the absence of boundary conditions for NDSolve[{{eq1,eq2},ics}//Flatten,...] you will notice that the plotted solution (last plot) does not fullfill at all v'(x)=0, as it would do for the second order telegraph equation with zero Neumann values. So I don't think that zero Neumann values are used in this precise case Mar 18, 2019 at 15:34