Skip to main content

differential equations - DSolve with SetOptions[Solve, Method -> Reduce] Behaving Differently in Ver 11.1.1


With version "11.1.1 for Microsoft Windows (64-bit) (April 18, 2017)",


DSolve[{SI'[s] == Sin[PH[s]] , PH'[s] == Sin[SI[s]]/e}, {SI[s], PH[s]}, s]


produces the common warning message,



Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information.



and then a lengthy answer, the Last half of which is given in the answer to 147688, all in about 5 sec. A typical method of employing Reduce, as the warning suggests, is


SetOptions[Solve, Method -> Reduce];

However, DSolve appears to ignore this setting, returning the same warning message and answer in the same amount of time.


With version "10.4.1 for Microsoft Windows (64-bit) (April 11, 2016)", and the default Method for Solve, DSolve again produces the same warning message and answer in the same amount of time. In contrast, with Method -> Reduce DSolve runs about 250 sec before returning unevaluated.



Question - How is this difference to be interpreted? Possibilities that come to mind include:



  1. Version 11.1.1 of DSolve immediately recognizes that Reduce cannot produce an answer and ignores the instruction to use it.

  2. Version 11.1.1 contains a bug, albeit fortuitous in this case.

  3. There is a new way to invoke Reduce with DSolve, and I am the last to know it.


By the way, this is not the first instance in which I have seen DSolve ver 11.1.1 ignore Method -> Reduce.



Answer



Introduction


I will discuss the issues that lead to the inverse function message and a few debugging tools I used, which it might be well to have collected in one place for reference. Finally I will show how to use some of these tools on the OP's example. In getting around the inverse function warning we are able to recover a branch winding number as a discrete parameter that is missed by the straight DSolve call. With respect to the main point, the question itself about how to control the method used by Solve is addressed. In fact how to target the specific call to Solve that emits the Solve::ifun message is shown.



Another thing to keep in mind is that DSolve analyzes the differential equation for particular special types. It is conceivable that a message might be generated in checking a submethod that is irrelevant (i.e. that would not solve it). Fixing the error message would not help in such a case.


Discussion of the problem


One can see (using techniques discussed in the following sections) that the offending equation seems quite simple:


Solve[-Cos[SI[PH]]==C[1]-e Cos[PH],SI[PH]]

If you try to solve it outright, there is no problem. You get the following with no error/warning messages:


{{SI[PH] ->
ConditionalExpression[-ArcCos[-C[1] + e Cos[PH]] + 2 π C[2], C[2] ∈ Integers]},
{SI[PH] ->
ConditionalExpression[ArcCos[-C[1] + e Cos[PH]] + 2 π C[2], C[2] ∈ Integers]}}


Note the term 2 π C[2], where the integer parameter C[2] could be called a "winding number." It turns out that DSolve[] does two things that affect the ability of Solve[] to handle the equation. It makes the following calls:


SetSystemOptions[ReduceOptions->UseTranscendentalSolve->False]
SetOptions[Solve, Method -> "Restricted"]

Note that the second one, SetOptions[], is executed after DSolve starts and overrides the SetOptions[Solve, Method -> Reduced] in this answer. On the other hand, it is "ReduceOptions" setting that leads to the Solve::ifun warning in the present case, which you can check with the following, with or without the Method setting:


With[{redopts = SystemOptions["ReduceOptions"]},
Internal`WithLocalSettings[
SetSystemOptions["ReduceOptions" -> "UseTranscendentalSolve" -> False],
Solve[-Cos[SI[PH]] == C[1] - e Cos[PH], SI[PH](*, Method -> "Restricted"*)],

SetSystemOptions[redopts]
]]

One thing to keep in mind is that automatic settings depend on heuristics, which are chosen to be helpful in most (but not always all) cases. I suspect that solutions containing ConditionalExpression are more trouble than they are worth, especially ones depending on a new parameter like C[2]. Hence, the "ReduceOptions" setting used by DSolve leads to the solution in the example. Likewise, the undocumented Method -> "Restricted", whatever it does, must remove or tweak cases that arise in the Automatic setting and are troublesome for DSolve.


Some debugging tools


It may be of general interest and useful to collect some of the ways of investigating DSolve. A few of these are undocumented but often found on this site. Some of these techniques can be used with other solvers, too. It seems almost obvious that a general solver DSolve might call other solvers like Integrate and Solve, and possibly others like Reduce (no) and Simplify (yes). It is worth remarking that while it appears from the docs that DSolve does not use assumptions, some of these other functions do, Integrate and Simplify being noteworthy examples. However, it does reset $Assumptions before returning (as I discovered below).


Tracing. There are Trace and On. On[Solve] and On[SetOptions] yields information on calls to each function. On[SetOptions] generates quite a few calls. A more targeted approach is to use Trace with patterns:


Trace[
DSolve[{SI'[s] == Sin[PH[s]], PH'[s] == Sin[SI[s]]/e}, {SI[s], PH[s]}, s],
_Solve | SetOptions[Solve, _],

TraceInternal -> True] // Flatten

Trapping system function calls. Internal`InheritedBlock allows one to intercept system function calls, add definitions to them and print out or Sow information about the call. The great thing here is that when it is done, the system function reverts to its old self, and you don't have to worry about forgetting to reset something.


Villegas-Gayley. The Villegas-Gayley trick is the way to insert your own code between a system function call and the execution of the function.


f[args...] /; !TrueQ[someVariable] := Block[{someVariable = True},
;
f[args]]

You have to get the new definition ahead of all others, and sometimes one has to do it explicitly by setting its DownValues.


Temporarily setting system options (or whatever). A typical Internal`WithLocalSettings usage has the form



With[{currentsettings = ...},
Internal`WithLocalSettings[
,
,

]]

The is executed first, then , then . The advantage here is that the reset happens even if the is aborted (e.g., by the menu command Evaluation > Abort Evaluation).


Intercepting Solve[] and solving the differential equation


This section serves two functions, to compute an enhanced solution and to show how to alter the options of Solve when called by DSolve. We will use many of the things discussed above to reset SystemOptions and change the options passed to Solve. In particular, Check is used to make an alternative call to Solve when the original call emits the Solve:ifun message.



In fact the alternate call to Solve returns conditional expressions, which DSolve did not handle successfully nor quickly. So I thought to break it down and add the condition to $Assumptions and post-proces the result of DSolve to be a ConditionalExpression based on the conditions accumulated. More importantly in this case, the conditions depend on a integer parameter that with the default settings is C[2], the same as a continuous parameter used by DSolve in the general solution. This conflict has to be resolved, which I did by manually setting the GeneratedParameters option in Solve. Another potential bug is that independent calls to Solve might generate independent parameters with the same name. While it in fact happens in this case, ignoring it turns out to be safe. (It could be dealt with, but this example does not seem worth the trouble.)


A utility to add conditions to a DSolve solution:


(* convert a DSolve solution to a ConditionExpression *)
addConditions[conditions_] := {
HoldPattern[var_ -> Function[x_, expr_]] :>
var -> With[{cond = Simplify[And @@ Flatten@conditions]},
Function @@ Hold[x, ConditionalExpression[expr, cond]]],
HoldPattern[var_ -> expr_] :>
var -> ConditionalExpression[expr, Simplify[And @@ Flatten@conditions]]};


The example intercept of Solve within DSolve:


On[$Assumptions];  (* shows when  $Assumptions  is changed *)
Module[{assum = True},
Internal`InheritedBlock[{Solve, $Assumptions},
Unprotect[Solve];
call : Solve[eq_, v_, opts___] /; ! TrueQ[$in] :=
Block[{$in = True, $res1, $res2},
Check[ (* Use Check to try default settings & respond to messages *)
$res1 = call (* try original call *)
, (* Check[] *)

Print["Trying \"ReduceOptions\" -> \"UseTranscendentalSolve\" -> True on ",
HoldForm[call]];
With[{redopts = SystemOptions["ReduceOptions"]},
Internal`WithLocalSettings[
SetSystemOptions[
"ReduceOptions" -> "UseTranscendentalSolve" -> True], (* I`WLS: init *)
$res2 = Solve[eq, v, GeneratedParameters -> $S, opts];
$res2 = $res2 /. ConditionalExpression[e_, c_] :> (* map condition to assumption *)
($Assumptions = $Assumptions && c; (* for solvers and Simplify[] *)
assum = assum && c; (* for addConditions[] at end *)

e);
Print["$Assumptions now ", $Assumptions]; (* context for On[$Assumptions] *)
$res2 = Simplify[$res2];
Print["Result= ", $res2],
SetSystemOptions[redopts] (* I`WLS: reset *)
];
];
If[FreeQ[$res2, Solve], $res2, $res1]
, (* Check[] *)
{Solve::ifun}]

];
Protect[Solve];
sol = DSolve[{SI'[s] == Sin[PH[s]], PH'[s] == Sin[SI[s]]/e}, {SI,
PH}, s];
Print["$Assumptions finally ", $Assumptions]; (* context for On[$Assumptions] *)
sol = sol /. addConditions[assum]
]]
Off[$Assumptions];

Message and Print output (Print statements shown in italics):




$Assumptions::trace: $Assumptions --> True. ... (x3)


Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information.


Trying "ReduceOptions" -> "UseTranscendentalSolve"->True on Solve[-Cos[SI[PH]]==C[1]-e Cos[PH],SI[PH]]


$Assumptions::trace: $Assumptions --> True.


$Assumptions::trace: $Assumptions --> $S[1]∈Integers.


$Assumptions::trace: $Assumptions --> $S[1]∈Integers && $S[1]∈Integers.


$Assumptions now $S[1]∈Integers && $S[1]∈Integers


$Assumptions::trace: $Assumptions --> $S[1]∈Integers && $S[1]∈Integers.


Result= {{SI[PH]->-ArcCos[-C[1]+e Cos[PH]]+2 π $S[1]},{SI[PH]->ArcCos[-C[1]+e Cos[PH]]+2 π $S[1]}}



$Assumptions::trace: $Assumptions --> $S[1]∈Integers && $S[1]∈Integers.


$Assumptions::trace: $Assumptions --> True.


$Assumptions finally True


$Assumptions::trace: $Assumptions --> True. ... (x4)



Output (abridged). The $S[1] is the "winding number" parameter lost in the default DSolve by "Inverse functions...being used"; it appears only in SI.


{{SI -> Function[{s},
ConditionalExpression[-ArcCos[ ..] + 2 π $S[1], $S[1] ∈ Integers]],
PH -> Function[{s},
ConditionalExpression[InverseFunction[.. &][-(s/e) + C[2]], $S[1] ∈ Integers]]},

{SI -> Function[{s},
ConditionalExpression[ArcCos[ ..] + 2 π $S[1], $S[1] ∈ Integers]],
PH -> Function[{s},
ConditionalExpression[InverseFunction[.. &][s/e + C[2]], $S[1] ∈ Integers]]}}

Example print statements and Solve commands that might be used to investigate :


(* Print statements *)
Print[HoldForm[call]];
Print["SystemOptions[\"ReduceOptions\"]= ", SystemOptions["ReduceOptions"]];
Print["opts= ", {opts}];

Print["Trying Method -> Reduce on ", HoldForm[call],
", Options[Solve]= ", Options[Solve],
", SystemOptions[\"ReduceOptions\"]= ",
SystemOptions["ReduceOptions"]];
Print["Trying \"ReduceOptions\" -> \"UseTranscendentalSolve\" -> True on ",
HoldForm[call]];
Print["Result= ", $res2];
Print@Stack[];

(* Alternative Solve[] ideas *)

Solve[eq, v, GeneratedParameters -> $S, Method -> Reduce, opts] (* Method -> Reduce *)
TimeConstrained[Solve[eq, v, GeneratedParameters -> $S, opts], 10] (* abort long Solve *)

References for undocumented features


Villegas-Gayley: What are some advanced uses for Block?, SO/4198961.


Internal`InheritedBlock: How safe is the use of Block and Internal`InheritedBlock


Internal`WithLocalSettings: SO/7586735, How to flush machine underflows to zero and prevent conversion to arbitrary precision?


Comments

Popular posts from this blog

plotting - Filling between two spheres in SphericalPlot3D

Manipulate[ SphericalPlot3D[{1, 2 - n}, {θ, 0, Pi}, {ϕ, 0, 1.5 Pi}, Mesh -> None, PlotPoints -> 15, PlotRange -> {-2.2, 2.2}], {n, 0, 1}] I cant' seem to be able to make a filling between two spheres. I've already tried the obvious Filling -> {1 -> {2}} but Mathematica doesn't seem to like that option. Is there any easy way around this or ... Answer There is no built-in filling in SphericalPlot3D . One option is to use ParametricPlot3D to draw the surfaces between the two shells: Manipulate[ Show[SphericalPlot3D[{1, 2 - n}, {θ, 0, Pi}, {ϕ, 0, 1.5 Pi}, PlotPoints -> 15, PlotRange -> {-2.2, 2.2}], ParametricPlot3D[{ r {Sin[t] Cos[1.5 Pi], Sin[t] Sin[1.5 Pi], Cos[t]}, r {Sin[t] Cos[0 Pi], Sin[t] Sin[0 Pi], Cos[t]}}, {r, 1, 2 - n}, {t, 0, Pi}, PlotStyle -> Yellow, Mesh -> {2, 15}]], {n, 0, 1}]

plotting - Plot 4D data with color as 4th dimension

I have a list of 4D data (x position, y position, amplitude, wavelength). I want to plot x, y, and amplitude on a 3D plot and have the color of the points correspond to the wavelength. I have seen many examples using functions to define color but my wavelength cannot be expressed by an analytic function. Is there a simple way to do this? Answer Here a another possible way to visualize 4D data: data = Flatten[Table[{x, y, x^2 + y^2, Sin[x - y]}, {x, -Pi, Pi,Pi/10}, {y,-Pi,Pi, Pi/10}], 1]; You can use the function Point along with VertexColors . Now the points are places using the first three elements and the color is determined by the fourth. In this case I used Hue, but you can use whatever you prefer. Graphics3D[ Point[data[[All, 1 ;; 3]], VertexColors -> Hue /@ data[[All, 4]]], Axes -> True, BoxRatios -> {1, 1, 1/GoldenRatio}]

plotting - Mathematica: 3D plot based on combined 2D graphs

I have several sigmoidal fits to 3 different datasets, with mean fit predictions plus the 95% confidence limits (not symmetrical around the mean) and the actual data. I would now like to show these different 2D plots projected in 3D as in but then using proper perspective. In the link here they give some solutions to combine the plots using isometric perspective, but I would like to use proper 3 point perspective. Any thoughts? Also any way to show the mean points per time point for each series plus or minus the standard error on the mean would be cool too, either using points+vertical bars, or using spheres plus tubes. Below are some test data and the fit function I am using. Note that I am working on a logit(proportion) scale and that the final vertical scale is Log10(percentage). (* some test data *) data = Table[Null, {i, 4}]; data[[1]] = {{1, -5.8}, {2, -5.4}, {3, -0.8}, {4, -0.2}, {5, 4.6}, {1, -6.4}, {2, -5.6}, {3, -0.7}, {4, 0.04}, {5, 1.0}, {1, -6.8}, {2, -4.7}, {3, -1.