# Obtain a stable solution using NDSolve

I am trying to solve a partial differential equation using NDSolve. When I am not specifying any method options with NDSolve, it gives me a solution which is quite oke but shows oscilations that are not suppose to be there. I expected to be able to solve this using options like MaxStepSize to get rid off or at least reduce this oscilations. However, unfortunately, when I use for example MaxStepSize -> 0.01 I get results that do not make sense and errors like NDSolve::mconly: For the method IDA, only machine real code is available. Unable to continue with complex values or beyond floating-point exceptions. >> I do not really understand why I get this error. Hopefully someone could illuminate this to me. Also I guess MaxStepSize is not the way to give me a better solution, are there any other options to improve my solution or should I just be satisfied with the one I have so far?

Simplification of the code:

tmin = 0; tmax = 70;
xmin = 0; xmax = 183;
rate = 0.00084;
a = 96;
b = 5/3;

f[t_] :=  If[t > 1 && t < 30, rate, 0];

eq1 = i[x, t] == a *h[x, t]^b;
eq2 = D[i[x, t], x] + D[h[x, t], t] ==  f[t];
ic1 := h[x, 0] == 0;
ic2 := i[x, 0] == 0;
bc1 := h[0, t] == 0;

solution = NDSolve[{eq1, eq2, ic1, ic2, bc1},
{i, h}, {x, xmin, xmax}, {t, tmin, tmax}, (*MaxStepSize -> 0.01*)];
Manipulate[
Plot[Evaluate[i[x, t] /. solution], {t, 0, tmax}, PlotStyle -> Red,
ImageSize -> {350}], {x, xmin, xmax}]


-
Try changing f to f2[t_] = UnitStep[t - 1] UnitStep[30 - t] rate and removing the MaxStepSize option. –  b.gatessucks Sep 20 '13 at 9:58
Thanks for your response, however, I think it does not really change my results? I see now I put in MaxStepSize, in my original code, I did not, I just tried this option to reduce the osciliation –  Wiebe Sep 20 '13 at 10:01
I get no error messages and I don't see the oscillations (unless I misunderstood). –  b.gatessucks Sep 20 '13 at 10:02
Maybe I was not completely clear. When I put in the code, I also do not get the error message. I get (for example at x=16.2)the picture that is shown above. This picture is almost what I expect, except that the top is suppose to be a straight line and not fluctuating. So I thought this was caused due to nummerical errors, which I could perhaps solve using MaxStepSize? Though when decreasing this stepsize I do get errormessages. –  Wiebe Sep 20 '13 at 10:50
@user9022 what makes you think it cant oscilate like that? –  george2079 Sep 20 '13 at 17:57

I think the reason why you get that error message is the term h[x, t]^b with b=5/3. If h[x,t] becomes negative for whatever reason i[x,t] will become complex and thus your error message. My guess is that i[x,t] becomes negative for numeric errors/instability here.

As you are solving a partial differential equation you need to also decrease the spatial grid, you can't just reduce the time step size only, these must fit together otherwise you'll produce instabilities. Here is an example for settings (its from the documentation in "possible issues" -> "partial differential equations") that does reduce the oscillation for your example on my machine (64bit):

solution =
NDSolve[
{eq1, eq2, ic1, ic2, bc1}, {i, h}, {x, xmin, xmax}, {t, tmin, tmax},
MaxStepSize -> 0.5, AccuracyGoal -> 8, PrecisionGoal -> 8,
Method -> {"MethodOfLines","SpatialDiscretization" -> {
"TensorProductGrid", "MinPoints" -> 1000}
}
];


Of course that will take much more time and memory to solve. In many comparable situations a good idea is to understand how important (realistic?) the sharp cuts at t=1 and t=30 are as these naturally are difficult for numeric solvers to handle correctly. By using smoother conditions there you might be able to get a decent result with much less computing time and memory usage.

You should also understand that playing with these options might get you a correct or plausible result but a detailed understanding about whether your solution can be trusted might take some extra effort and can't be deduced just because Mathematica doesn't complain. As I didn't take that extra effort I can't say how reliable my settings really are (it is easy enough to create instabilities with just slight changes...).

When playing with these options it might also be a good idea to use something like MemoryConstrained[NDSolve[...],1000000] (adopt to the amount of memory of your machine) to keep your system reactive, at least on Windows...

-
Thank you very much for your clear explanation. Do you (or anyone else) have any advice on documentation I could read to learn how to check whether the solution can be trusted? In this case I knew the result on forehand, but if this corresponds to the model results, does this still not assure that the model is right? –  Wiebe Sep 23 '13 at 5:49
@user9022: investigation of the stability/accurateness of numerical algorithms for solving differential equations can become a difficult problem depending on the degree of rigor you are after. There are many textbooks about the topic and I'm not enough of an expert to make a good suggestion. You should probably also clarify what degree of rigor you need, for many practical cases you certainly can get away with plausability checks as "does the solution look OK". Whether the model is "right" (== are the equations a correct description of the "real system") is a completely different question... –  Albert Retey Sep 23 '13 at 9:51