Re: [Help-gsl] Re: (A newbie asks) How to program multidimensional root solver?

[John Pye]

Hi all

I realise that this is rather off-topic for GSL but perhaps it will be
of interest in any case.

Olumide wrote:
>>> Thanks for replying, no I haven't. But I will. I've just downloaded
>>> Maxima, and ASCEND and I'm taking a closer look at them. Which of
>>> these two programs do you think is best for solving system on
>>> non-linear equations? BTW, I don't expecting an analytic solution. My
>>> equations don't look too friendly ... (nasty logs! :-) )
>>>
>> If you just want a numerical solution then try ASCEND or GAMS (or AMPL
>> or...). The advantage here is that you can enter the equations in a
>> fairly intuitive way and not have to worry about 'wiring it up' to the
>> solver. The Excel Solver (or even maybe the new Open Office 3.0b
>> Solver?) would be a good bet too. Maxima is really for getting analytic
>> solutions, although it might be possible to get numerical solutions as
>> well -- I'm not sure.
>
> What do you think of scilab's fsolve? -- I'm assuming you're an
> expert-type person in this area ;-) . I'm not liking ASCEND's
> documentation. ('Way too much to fluff to wade though.)
>

I agree about the ASCEND documentation, it is hard to quickly get the
big picture of how it all works.

I have attached a model that will maybe get you started.

You asked about fsolve in Scilab. ASCEND has a lot of stuff to make
building up of large models much easier than it is with tools like
fsolve: stuff to do with the easier of stating the model in a concise
way, object-oriented mathematical models, dealing with units of
measurement, playing around with which variables you want to solve for,
etc. In terms of pure numerics, fsolve is fine, just as GSL is fine. I
should mention that I'm involved with the ASCEND project, so I'm a bit
biased :-)

See attached.

Cheers
JP

REQUIRE "atoms.a4l";

MODEL olumide;
        n IS_A integer_constant;
        n :== 2;

        a, b, c, d, e, f, g[0..n], h[0..n], K, w[0..n] IS_A solver_var;
        
        a + b*x + c*y + d*x^2 + e*x*y + f*y^2 + SUM[r[i]^2 * ln (r[i]*w[i]) | i 
IN [0..n]] = z;
        x + z*(b + 2*d*x + e*y + SUM[w[i]*(ln(r[i]) + K) *(x - g[i]) | i IN 
[0..n]]) = 0;
        y + z*(c + 2*f*y + e*x + SUM[w[i]*(ln(r[i]) + K)*(x - h[i]) | i IN 
[0..n]]) = 0;

        x, y, z, r[0..n] IS_A solver_var;
        
        FOR i IN [0..n] CREATE
                r[i] = ( (x - g[i])^2 + (y - h[0])^2 )^0.5;
        END FOR;

METHODS
METHOD on_load;
        FIX a, b, c, d, e, f, g[0..n], h[0..n], K, w[0..n];
        a := 1;
        b := 1;
        c := 1;
        d := 1;
        e := 1;
        f := 1;
        g[0] := 1;
        g[1] := 4;
        g[2] := 1;
        h[0] := 1;
        h[1] := -4;
        h[2] := 4;
        K := 1;
        w[0] := 1;
        w[1] := 1;
        w[2] := 1;

END on_load;
END olumide;