2
$\begingroup$

The following piece of code

     (* DEFINITIONS *)
    mol[n_Integer, o_: "Pseudospectral"] := {"MethodOfLines",
      "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n, 
        "MinPoints" -> n, "DifferenceOrder" -> o}}  


    f[t_] := Exp[t]/(1 + Exp[t])
    step[t_] := f[(t - 1/2)*(100)] - f[(0 - 1/2)*(100)]
    (*  I tried to make sure that step[t] is exactly zero at t=0 to avoid initial violation of the BC's*)



    (* PARAMETERS *)
    k = 1/3; 
    A = 1/40; 
    yo = 3/10;


    (* INTEGRATION*)
    s = NDSolve[{
        D[v[t, x], t] + v[t, x]*D[v[t, x], x] + k*D[y[t, x], x] ==
          0,
        D[y[t, x], t] + v[t, x]*D[y[t, x], x] +       D[v[t, x], x] == 0,

        y[0, x] == yo,
        y[t, 0] == yo,

        v[0, x] == 0,
        v[t, 1] == A*step[t]


        }, {v, y}, {x, 0, 1}, {t, 0, 1}, Method -> mol[81]];

    (* VISUALIZATION *)
    Plot[(v[t, 1] - A*step[t]) /. s, {t, 0, 1}, WorkingPrecision -> 100, 
     PlotPoints -> 400, PlotRange -> All, PlotLabel -> "BC violation"]

results in a plot that shows that BC are violated on their 3rd decimal digit.

Is there any typo or spelling mistake in my code? How can I fix it?

APPENDIX: Another example of violation that gets close to the second decimal digit

    (* DEFINITIONS *)
mol[n_Integer, o_: "Pseudospectral"] := {"MethodOfLines",
  "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n, 
    "MinPoints" -> n, "DifferenceOrder" -> o}}  


f[t_] := Exp[t]/(1 + Exp[t])
step[t_] := f[(t - 1/2)*(100)] - f[(0 - 1/2)*(100)]


(* PARAMETERS *)
λ = 1/3; 
A = 1/6; 
yo = 2/10;


(* INTEGRATION*)
s = NDSolve[{
    D[v[t, x], t] + v[t, x]*D[v[t, x], x] + λ*D[y[t, x], x] ==
      0,
    D[y[t, x], t] + v[t, x]*D[y[t, x], x] +       D[v[t, x], x] == 0,

    y[0, x] == yo,
    y[t, 0] == yo + A*step[t],

    v[0, x] == 0,
    v[t, 1] == 0


    }, {v, y}, {x, 0, 1}, {t, 0, 1}(*,Method\[Rule]mol[81]*),
   Method -> {"MethodOfLines",
     "SpatialDiscretization" -> {"TensorProductGrid", 
       "MaxPoints" -> 2001, "MinPoints" -> 2001, 
       "DifferenceOrder" -> Automatic}}  ];

(* VISUALIZATION *)
Plot[(y[t, 0] - yo - A*step[t]) /. s, {t, 0, 1}, 
 WorkingPrecision -> 100, PlotPoints -> 400, PlotRange -> All, 
 PlotLabel -> "BC violation"]
Plot[(v[t, 1]) /. s, {t, 0, 1}, WorkingPrecision -> 100, 
 PlotPoints -> 400, PlotRange -> All, PlotLabel -> "BC violation"]
$\endgroup$

1 Answer 1

1
$\begingroup$

This is related to how b.c. is imposed in NDSolve. In short, the b.c. has actually been rebuilt by numerically solving ODEs, thus numeric error comes up. A detailed discussion can be found here. To resolve the issue, you can make the step size in $t$ direction smaller:

s2 = NDSolve[{D[v[t, x], t] + v[t, x]*D[v[t, x], x] + k*D[y[t, x], x] == 0, 
    D[y[t, x], t] + v[t, x]*D[y[t, x], x] + D[v[t, x], x] == 0, y[0, x] == yo, 
    y[t, 0] == yo, v[0, x] == 0, v[t, 1] == A*step[t]}, {v, y}, {x, 0, 1}, {t, 0, 1}, 
   Method -> Union[mol[81]], MaxStepSize -> {0.001, Automatic}];


Plot[(v[t, 1] - A*step[t]) /. s2 // Evaluate, {t, 0, 1}, PlotRange -> All, 
 PlotLabel -> "BC violation"]

enter image description here

Turn to the DAE solver will make the result even better (but do notice DAE solver is generally weaker than the ODE solver, at least now):

mol[tf:False|True,sf_:Automatic]:={"MethodOfLines",
"DifferentiateBoundaryConditions"->{tf,"ScaleFactor"->sf}}

s3 = NDSolve[{D[v[t, x], t] + v[t, x]*D[v[t, x], x] + k*D[y[t, x], x] == 0, 
    D[y[t, x], t] + v[t, x]*D[y[t, x], x] + D[v[t, x], x] == 0, y[0, x] == yo, 
    y[t, 0] == yo, v[0, x] == 0, v[t, 1] == A*step[t]}, {v, y}, {x, 0, 1}, {t, 0, 1}, 
   Method -> Union[mol[81], mol[False]], MaxStepSize -> {0.001, Automatic}];

Plot[(v[t, 1] - A*step[t]) /. s3 // Evaluate, {t, 0, 1}, PlotRange -> All, 
 PlotLabel -> "BC violation"]

enter image description here

BTW, if you want to check the solution with Plot3D, remember to enlarge PlotPoints:

Plot3D[(v[t, x]) /. s3 // Evaluate, {t, 0, 1}, {x, 0, 1}, PlotRange -> All, 
 PlotPoints -> 50]

enter image description here

$\endgroup$

Your Answer

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