0

Consider the following numerical problem. One deals with transcendental complicated function f that depends on several parameters and three variables, say x, y and z. One wants to find roots of the equation f[x,y,z]==0 for a given values x1 and y1 (i.e. solve f[x,y,z] with respect to z for known values x and y. Important moment that it is possible to have more thatn one root of f[x,y,z]==0 for x1 and y1.

Then, for a given lists {x} and {y}, one solves equation at each point on grid (x,y) and obtain the list z. It can be realized as follows (schematically!):

roots = z/.NDSolve[f[#1,#2,z]==0 && g[#1,#2]<=z<=h[#1,#2]]&;
zvalues = ParallelTable[root[X,Y], ... ]

I write down ParallelTable because I am not sure that the function roots is listable, i.e. I am not sure that Map or MapThread works.

Then, one notices that this procedure is quite long. The first guess is to try use Compile.

So, the following questions appear:

  1. Is it possible to find numerically all roots of equation on a given interval by a method, that can be used in Compile?
  2. Or, is it possible to use NSolve in Compile?

I have a vague feeling that this question is relevant. With help of kind people, I obtain the following code snippet:

NewtonMethodCompiled::usage = 
"NewtonMethodCompiled[function][x0]. \n\
NewtonMethodCompiled[function, delta, epsilon][x0]. \n\
NewtonMethodCompiled[function, \"Delta\" -> delta, \"Epsilon\" -> epsilon][x0]. ";

Options[NewtonMethodCompiled] = {"Delta" -> 0.0001, "Epsilon" -> 0.0001};


NewtonMethodCompiled[f_Function, OptionsPattern[]] := 
NewtonMethodCompiled[f, OptionValue["Delta"], OptionValue["Epsilon"]];


NewtonMethodCompiled[f_Function, delta_Real, epsilon_Real] := 
NewtonMethodCompiled[f, delta, epsilon] = 
Compile[{{x0, _Complex}}, 
    Module[{xi, df}, 
        xi = x0; 
        While[epsilon < Abs[f[xi]] < Abs[f[x0]] / epsilon, 
            df = (f[xi + delta/2] - f[xi - delta/2]) / delta;
            xi = xi - f[xi] / df; 
        ];


    Return[xi] 
], 
CompilationOptions -&gt; {&quot;InlineExternalDefinitions&quot; -&gt; True, &quot;InlineCompiledFunctions&quot; -&gt; True}



];

It is the realization of compiled Newton method

Artem Alexandrov
  • 803
  • 4
  • 11
  • Compile can make external calls, so, technically yes, you can use such tools as NSolve, NDSolve, etc. as long as your code returns a type of expression supported by the WVM (machine numeric/boolean scalars and a arrays). It is possible to wrap NSolve in code that returns a numeric array of solutions. However, functions that don't return supported data types cannot be run solely within the WVM. -- I suppose you've probably seen this: https://mathematica.stackexchange.com/questions/1096/list-of-compilable-functions – Michael E2 Apr 29 '21 at 23:50
  • @MichaelE2 thank you! My function that contains NSolve returns list of reals and naive compilation produces CompiledFunction::cfse: Compiled expression {0.396782,0.785873} should be a machine-size real number. . I will check documentation and the mentioned question – Artem Alexandrov Apr 30 '21 at 06:21

0 Answers0