I have a problem with obtaining the domain of resulting Interpolating Function in my calculations. I want to solve the following differential equation

V[u_, v_] =ProductLog[E^(1/2 (-2 - u + v))]/(1 + ProductLog[E^(1/2 (-2 - u + v))])^4;
sol = First[NDSolve[{-4*D[S[u, v], u, v] == V[u, v]* S[u, v], 
S[u, 0] == Exp[-(100/18)], S[0, v] == Exp[-((v - 10)^2/18)]},S, {u,100,200}, {v,100,200},
Method -> {"MethodOfLines","SpatialDiscretization" -> { "TensorProductGrid", 
"MaxStepSize" -> 1}}, AccuracyGoal -> 1]]

As you can see, I determined the range of u and v from 100 to 200, but if you run the code, you will find that the resulting Interpolating Function range will be from 0 to 200 for v. How can I fix this problem? Thanks all.

  • $\begingroup$ It seems to be the presence of the statement S[u, 0] == Exp[-(100/18)] in your NDSolve command that places the lower bound on v. If you change that to S[u, 10] == ..., then the interpolating function changes accordingly. I can't speak to whether this is the intended behavior - it seems to treat the bounds on u and v differently. $\endgroup$ – Jason B. Aug 13 '18 at 20:48
  • $\begingroup$ @JasonB. Actually, S[u, v] == Exp[-(v - 10)^2/18] is the boundary condition of the problem which I used S[u, 0] and S[0, v] in calculations. Can I use S[u,100] == Exp[-(90)^2/18] as a new boundary condition to obtain the result that I am looking for? $\endgroup$ – Mehrab Aug 13 '18 at 21:16
  • $\begingroup$ Can't you just try the new BCs and see what happens? (I'm pretty sure that the reason for the v domain is that in the method of lines, the spatial variable's derivatives are approximated on a spatial grid between the boundaries. S[u, 0] specifies one boundary to be v == 0. Now why can't NDSolve discard all the points outside the specified interval? I don't know. I suppose one cannot have an arbitrary boundary coincide with a grid point. Maybe that's why.) $\endgroup$ – Michael E2 Aug 13 '18 at 23:01
  • 1
    $\begingroup$ @Mehrab Because the boundary conditions are at u == 0 and v == 0, NDSolve must integrate to those surfaces, whether or not it includes them in the InterpolatingFunction. That being the case, simply ignore the unwanted part of the InterpolatingFunction. By the way, "TemporalVariable" -> v yields a domain of {{0., 200.}, {100., 200.}}, and "TemporalVariable" -> u yields a domain of {{100., 200.}, {0., 200.}}. But, why worry about this idiosyncrasy? It does not seem to cause any harm. Moreover, it is easy enough to trim the domain by postprocessing, if desired. $\endgroup$ – bbgodfrey Aug 14 '18 at 3:43
  • 1
    $\begingroup$ @Mehrab I should add that AccuracyGoal -> 1 does not yield a very accurate answer. Just omit this option. $\endgroup$ – bbgodfrey Aug 14 '18 at 3:57

You can consider solving the problem in two steps, and then compare for different options AccuracyGoal. Below is the code and comparison for the four options AccuracyGoal. In Fig. it is clear that the oscillations of the solution are a consequence of the numerical error. These oscillations disappear as the accuracy increases.

V[u_, v_] = 
    E^(1/2 (-2 - u + v))]/(1 + ProductLog[E^(1/2 (-2 - u + v))])^4;
sol = First[
  NDSolve[{-4*D[S[u, v], u, v] == V[u, v]*S[u, v], 
    S[u, 0] == Exp[-(100/18)], S[0, v] == Exp[-((v - 10)^2/18)]}, 
   S, {u, 0, 200}, {v, 0, 200}, 
   Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"TensorProductGrid", 
       "MaxStepSize" -> 1}}, AccuracyGoal -> 15]]
S0[u_, v_] := S[u, v] /. sol
sol1 = First[
  NDSolve[{-4*D[S[u, v], u, v] == V[u, v]*S[u, v], 
    S[u, 100] == S0[u, 100], S[100, v] == S0[100, v]}, 
   S, {u, 100, 200}, {v, 100, 200}, 
   Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"TensorProductGrid", 
       "MaxStepSize" -> 1}}, AccuracyGoal -> 15]]
{Plot3D[Evaluate[S[u, v] /. sol1], {u, 100, 200}, {v, 100, 200}, 
  Mesh -> None, PlotRange -> All, PlotLabel -> ""], 
 Plot3D[Evaluate[S[u, v] /. sol], {u, 100, 200}, {v, 100, 200}, 
  Mesh -> None, PlotRange -> All]}


  • $\begingroup$ Thanks a lot. Although you answered the question, unfortunately, my problem is still there. For more details please have a look at my comments to Godfrey given above. $\endgroup$ – Mehrab Aug 14 '18 at 14:50
  • $\begingroup$ Mehrab, you can not break the 2D region into intervals, you can break the region into subdomains. Thus, there is no benefit of partitioning into intervals, since it is required to specify a function on the boundaries of the domain, and not on the boundaries of the interval. $\endgroup$ – Alex Trounev Aug 14 '18 at 15:08
  • $\begingroup$ Ok, thanks. I guess I should look for a supercomputer for my calculation and there is no other choice. $\endgroup$ – Mehrab Aug 14 '18 at 15:31
  • $\begingroup$ A supercomputer is not enough to solve a problem. We need a good algorithm. With a good algorithm, we do not need a supercomputer. I'll see what else can be done in this task. $\endgroup$ – Alex Trounev Aug 14 '18 at 17:17
  • $\begingroup$ Yes, you are right. Thanks for your kind effort. How can I become an expert in Mathematica? Do I need to pass some courses? Or I just need to read a book. Do you have any suggestion? $\endgroup$ – Mehrab Aug 14 '18 at 18:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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