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 = {#, 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}]]


v[0] in dashed, v[t] in blue, i[t] in cyan

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]];

v[0] in dashed, v[t] in blue, i[t] in cyan

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]];

unspecified  boundary conditions

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.

  • $\begingroup$ Do you want to find a solution for non-reflecting boundary conditions? $\endgroup$ Mar 14, 2019 at 15:40
  • $\begingroup$ 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. $\endgroup$ Mar 15, 2019 at 9:35

1 Answer 1


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]


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}]

enter image description here

 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.

  • $\begingroup$ Thank you user 21 $\endgroup$ Mar 15, 2019 at 23:03
  • $\begingroup$ 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. $\endgroup$ Mar 15, 2019 at 23:13
  • $\begingroup$ 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}? $\endgroup$ Mar 18, 2019 at 15:26
  • $\begingroup$ And actually, I want non-reflective conditions on both ends (no Dirichlet condition) but it is a detail. $\endgroup$ Mar 18, 2019 at 15:28
  • $\begingroup$ 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 $\endgroup$ Mar 18, 2019 at 15:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.