We have 2 (3) ways of proceed :
make experiments, find good
expert on numerical methods (do both)
The script attached shows, for small x
the best h is as big as 10^5*sqrt($macheps)
Still I do not know what to to with x: abs(x)>10^5
and what to do with say, exp(100)
try
x = {100}
eval fdjac (x,exp(x)) - exp(100)
That's ERROR!
Oleh


30 грудня 2015, 03:08:30, від "Allin Cottrell" <cottrell@wfu.edu>:

On Tue, 29 Dec 2015, Riccardo (Jack) Lucchetti wrote:

> On Tue, 29 Dec 2015, Allin Cottrell wrote:
>
>> In current git I've modified fdjac2 to use a step size as in your script 
>> (when Richardson is selected). We need to do more testing on this in case 
>> it goes nasty under some conditions
>
> It would seem that the results are somewhat sensitive to the amount of 
> nonlinearity of the problem: I ran an old script I had which uses matrix 
> inversion (notoriously differentiable almost everywhere but highly nonlinear 
> --- especially in the nieighbourhood of singularity).
>
> <hansl>
> set echo off
> set messages off
>
> function matrix vec22inv(matrix x)
>    n = round(sqrt(rows(x)))
>    A = mshape(x,n,n)
>    return vec(inv(A))'
> end function
>
> function matrix grad(matrix x)
>    n = round(sqrt(rows(x)))
>    iA = inv(mshape(x,n,n))
>    return -iA' ** iA
> end function
>
> n = 12
> k = 1.0e4
> a = vec(ones(n,n) - I(n)/k)
>
> ia = vec22inv(a)
> ag = grad(a)
>
> printf "\n"
>
> loop i=0..2 --quiet
>    set fdjac_quality $i
>    set stopwatch
>    ng$i = fdjac(a, "vec22inv(a)")
>    tt = $stopwatch
>    scalar aae = meanc(abs(vec(ag - ng$i)))
>    scalar mae = maxc(abs(vec(ag - ng$i)))
>    scalar are = meanc(abs(1 - vec(ag ./ ng$i)))
>    scalar mre = maxc(abs(1 - vec(ag ./ ng$i)))
>    printf "quality = $i\n"
>    printf "\tAverage and max absolute error = %g, %g\n", aae, mae
>    printf "\tAverage and max relative error = %g, %g\n", are, mre
>    printf "\t(cpu time = %g)\n", tt
> endloop
> </hansl>
>
> This is what I get before and after the patch:
>
> <before>
> quality = 0
> 	Average and max absolute error = 53.1382, 11477.9
> 	Average and max relative error = 2.27644e-05, 0.00013658
> 	(cpu time = 0.02)
> quality = 1
> 	Average and max absolute error = 0.0163915, 6.19778
> 	Average and max relative error = 7.58353e-09, 1.68055e-07
> 	(cpu time = 0.04)
> quality = 2
> 	Average and max absolute error = 0.00128204, 0.149397
> 	Average and max relative error = 1.11463e-09, 2.15128e-07
> 	(cpu time = 0.09)
> </before>
>
> <after>
> quality = 0
> 	Average and max absolute error = 53.1382, 11477.9
> 	Average and max relative error = 2.27644e-05, 0.00013658
> 	(cpu time = 0.02)
> quality = 1
> 	Average and max absolute error = 0.0163915, 6.19778
> 	Average and max relative error = 7.58353e-09, 1.68055e-07
> 	(cpu time = 0.08)
> quality = 2
> 	Average and max absolute error = 298.462, 128846
> 	Average and max relative error = 0.000128065, 0.00153573
> 	(cpu time = 0.07)
> </after>

Interesting. The "after" case is pretty compelling evidence that 
Oleh's setting of the step size -- which worked so nicely in some 
cases -- can produce very bad results in others.

However, here (building with gcc 5.3 on Fedora 23), the results are 
not so impressive for the "quality 2" case under the status quo ante 
(git commit ee7808 of fdjac2.c). I'm seeing:

<before>
quality = 0
 	Average and max absolute error = 53.1431, 11477.9
 	Average and max relative error = 2.27665e-05, 0.000136579
 	(cpu time = 0.0192254)
quality = 1
 	Average and max absolute error = 0.0448263, 6.1979
 	Average and max relative error = 3.92204e-08, 8.20507e-07
 	(cpu time = 0.0307975)
quality = 2
 	Average and max absolute error = 0.0814845, 1.41264
 	Average and max relative error = 9.3409e-08, 2.03417e-06
 	(cpu time = 0.073777)
</before>

Do we have a case where "quality 2" unambiguously produces better 
results than "quality 1"?

Allin


_______________________________________________
Gretl-devel mailing list
Gretl-devel@lists.wfu.edu
http://lists.wfu.edu/mailman/listinfo/gretl-devel