[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: fsolve regression
From: |
Jaroslav Hajek |
Subject: |
Re: fsolve regression |
Date: |
Tue, 24 Feb 2009 14:04:34 +0100 |
On Fri, Feb 20, 2009 at 7:36 PM, Jaroslav Hajek <address@hidden> wrote:
> On Fri, Feb 20, 2009 at 7:18 PM, John W. Eaton <address@hidden> wrote:
>> The following script shows a problem with the new fsolve. It doesn't
>> seem to terminate for me (I interrupted it after several minutes).
>> After looking closer at the values passed to the "inflate" function,
>> it appears to be cycling among a few input values.
>>
>> The initial guess and expected result are
>>
>> z0 =
>>
>> 0.133000000000000
>> 0.325115654084158
>>
>> z =
>>
>> 0.134806285287651
>> 0.323387263653996
>>
>>
>> so the initial guess is fairly good.
>>
>> jwe
>>
>>
>>
>> load foo.dat
>>
>> more off
>>
>> global K1 K2 P ext0 Gval dels unk zlast
>>
>> function retval = G (var)
>> global K1 K2 P Gval ext0 unk
>> w = ext0;
>> w(unk) = var;
>> x = w(1);
>> y = w(2);
>> z = 0.5 - x - y;
>> p = 0.5 + z;
>> if (x == 0)
>> xlogx = 0;
>> else
>> xlogx = abs(x)*log(abs(x/p));
>> endif
>> if (y == 0)
>> ylogy = 0;
>> else
>> ylogy = abs(y)*log(abs(y/p));
>> endif
>> if (z == 0)
>> zlogz = 0;
>> else
>> zlogz = abs(z)*log(abs(z/p));
>> endif
>> retval = -Gval -(x*log(K1) + y*log(K2)) + p*log(P) + 2*zlogz + xlogx +
>> ylogy;
>> endfunction
>>
>> function retval = inflate (z)
>> global zlast dels
>> retval = [G(z); norm(z-zlast) - dels^2];
>> endfunction
>>
>> [z, fval, info] = fsolve (@inflate, z0);
>>
>> format long
>> z0
>> z
>>
>>
>>
>
> OK, thanks for the example. I'll try to investigate the problem.
>
>
Sorry for the delay. It was a trivial problem after all caused by not
updating iteration count after unsuccesful iteration, which caused
trust-region radius to be reinitialized over and over.
I'm not sure, however, if this is compatible with Matlab - does it
count unsuccesful iterations? I hope it does not matter to anyone. If
Matlab uses a different enough algorithm, it may not even make sense.
With the new patch, your script gives:
z0 =
0.133000000000000
0.325115654084158
z =
0.134806289085259
0.323387260675138
You can get even closer by setting lower TolX and TolFun (by default
they're at sqrt(eps), ~eps requests maximum precision).
regards
--
RNDr. Jaroslav Hajek
computing expert
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz