[Date Index]
[Thread Index]
[Author Index]
Re: Can I solve this system of nonlinear equations?
*To*: mathgroup at smc.vnet.net
*Subject*: [mg125278] Re: Can I solve this system of nonlinear equations?
*From*: "Stephen Luttrell" <steve at _removemefirst_stephenluttrell.com>
*Date*: Sun, 4 Mar 2012 04:33:06 -0500 (EST)
*Delivered-to*: l-mathgroup@mail-archive0.wolfram.com
*References*: <jil5ld$rrm$1@smc.vnet.net> <jit112$e7i$1@smc.vnet.net>
"Ray Koopman" <koopman at sfu.ca> wrote in message
news:jit112$e7i$1 at smc.vnet.net...
> On Feb 29, 4:28 am, Andy <andy7... at gmail.com> wrote:
>> I'm dealing with systems of nonlinear equations that have 8 equations
>> and 8 unknowns. Here's an example:
>>
>> Solve[{(((c - a)/0.002) - (0.995018769272803 + h*b)) == 0,
>> (((d - b)/0.002) - (0.990074756047929 + h*c)) == 0,
>> (((e - c)/0.002) - (0.985167483257382 + h*d)) == 0,
>> (((f - d)/0.002) - (0.980296479563062 + h*e)) == 0,
>> (((g - e)/0.002) - (0.975461279165159 + h*f)) == 0,
>> (((-1*e + 8*d - 8*b + a)/(12*0.001)) - (0.990074756047929 + h*c)) ==
>> 0,
>> (((-1*f + 8*e - 8*c + b)/(12*0.001)) - (0.985167483257382 + h*d)) ==
>> 0,
>> (((-1*g + 8*f - 8*d + c)/(12*0.001)) - (0.980296479563062 + h*e)) ==
>> 0}, {a, b, c, d, e, f, g, h}]
>>
>> Whenever I try this, Mathematica 7 just returns the empty set {}. How
>> can I tell if this is unsolvable? Shouldn't I at least be able to get
>> a numerical approximation with NSolve? I've tried using stochastic
>> optimization to get approximate answers but every method gives poor
>> results, and that's why I would like to at least approximately solve
>> this if possible. Thanks very much for any help~
>
> To see if there is an exact solution to the problem,
> make all the coefficients exact and use Solve:
>
> x = { 500(c - a) - (995018769272803*^-15 + h*b),
> 500(d - b) - (990074756047929*^-15 + h*c),
> 500(e - c) - (985167483257382*^-15 + h*d),
> 500(f - d) - (980296479563062*^-15 + h*e),
> 500(g - e) - (975461279165159*^-15 + h*f),
> (-1*e + 8*d - 8*b + a)1000/12 - (990074756047929*^-15 + h*c),
> (-1*f + 8*e - 8*c + b)1000/12 - (985167483257382*^-15 + h*d),
> (-1*g + 8*f - 8*d + c)1000/12 - (980296479563062*^-15 + h*e)};
>
> Solve[Thread[x == 0],{a,b,c,d,e,f,g,h}]
>
> {}
>
> To get an approximate solution, minimize some measure of the
> differences between the two sides of the equations. The sum of
> the squared differences is convenient and not obviously improper;
> that is, there is nothing to suggest that the equations are on
> vastly different scales.
>
> TimeConstrained[Minimize[x.x,{a,b,c,d,e,f,g,h}],300] aborted.
>
> Both Minimize[N[x.x],{a,b,c,d,e,f,g,h}] and
> NMinimize[x.x,{a,b,c,d,e,f,g,h}] finished in < 2 sec
> and gave identical results:
>
> {3.7850320584543324`*^-11,
> {a -> 1.78322069825487`,
> b -> 1.7856155567256693`,
> c -> 1.7880252569409105`,
> d -> 1.79041402958568`,
> e -> 1.7928176848807527`,
> f -> 1.7952004996719013`,
> g -> 1.7975982365756977`,
> h -> 0.7881096531282696`}}
>
> The fit is a little better than Stephen Luttrell got,
> but the parameter estimates are far what he got.
>
> x /. %[[2]]
>
> {-2.832840262367853`*^-7,
> 1.7089251762580915`*^-6,
> 3.9068509285478115`*^-6,
> 1.6397939603951528`*^-6,
> -2.7478474384778906`*^-7,
> -1.6998965117753784`*^-6,
> -3.348694662008711`*^-6,
> -1.6487347667126784`*^-6}
>
> It looks like only the first 5 digits of the 15-digit constants
> are accurate.
>
When I try to reproduce your solution using NMinimize I get my solution (not
yours). I am using version 8.0.1 for Windows (32-bit). I presume that the
discrepency is because the solution space is not a single point, so
different versions and/or methods will arrive at different points in the
solution space.
Anyway, I derived the Hessian matrix to have a close look at the
neighbourhood of my solution point to see how the objective function
behaved, where the solution space seems to be 2-dimensional. This doesn't
tell us what it might look like further away. More work would be needed to
characterise it all.
Equations.
eqns = {(((c - a)/0.002) - (0.995018769272803 + h*b)) ==
0, (((d - b)/0.002) - (0.990074756047929 + h*c)) ==
0, (((e - c)/0.002) - (0.985167483257382 + h*d)) ==
0, (((f - d)/0.002) - (0.980296479563062 + h*e)) ==
0, (((g - e)/0.002) - (0.975461279165159 + h*f)) ==
0, (((-1*e + 8*d - 8*b + a)/(12*0.001)) - (0.990074756047929 +
h*c)) ==
0, (((-1*f + 8*e - 8*c + b)/(12*0.001)) - (0.985167483257382 +
h*d)) ==
0, (((-1*g + 8*f - 8*d + c)/(12*0.001)) - (0.980296479563062 +
h*e)) == 0};
Objective function.
obj = Total[eqns /. (x : _) == 0 :> x^2];
Minimise.
vars = {a, b, c, d, e, f, g, h};
soln = NMinimize[obj, vars]
giving what I got before ...
{3.94032*10^-11, {a -> 0.314882, b -> 0.315979, c -> 0.317109,
d -> 0.318197, e -> 0.319318, f -> 0.320396, g -> 0.321509,
h -> 0.374577}}
Now compute the matrix of second derivatives (Hessian) of the objective
function at this solution point.
hessian = Outer[D[obj, #1, #2] &, vars, vars] /. soln[[2]];
Then compute the eigenvalues of the Hessian.
evals = Eigenvalues[hessian]
{4.19009*10^6, 3.42241*10^6, 1.39642*10^6,
1.00435*10^6, 403393., 0.0159205, -2.57089*10^-10, 7.18138*10^-12}
Eyeballing this (or using ListPlot[evals, Filling -> Axis]) you see 5 large
(order 10^6), 1 very small (order 10^-2), and 2 noise level (order 10^(-10
or -11)) eigenvalues, so the the objective function stays at its minimum
value as you move away from the solution point in 2 (or 3, almost)
independent directions. The solution space is 2-dimensional (at least)
around the solution point.
You can use EigenSystem[hessian] to look at this 2D solution space in more
detail.
{evals, evects} = Eigensystem[hessian];
The last 2 eigenvectors (i.e. with noise level eigenvalues) span the
2-dimensional solution space using coefficients (u,v).
solnoffset = Thread[vars -> vars + {u, v}.Take[evects, -2]]
{a -> a + 0.37754 u - 0.000830755 v,
b -> b + 0.377681 u - 0.000462092 v,
c -> c + 0.377823 u - 0.000199144 v,
d -> d + 0.377964 u + 0.000171975 v,
e -> e + 0.378106 u + 0.000437376 v,
f -> f + 0.378247 u + 0.000810937 v,
g -> g + 0.378389 u + 0.00107878 v,
h -> h - 0.000381899 u + 0.999999 v}
Check what the objective function now looks like in the 2-dimensional
neighbourhood of the solution point.
objsoln[u_, v_] = obj /. solnoffset /. soln[[2]] // Expand // Chop
1.66682*10^-7 u^4 - 0.000872909 u^3 v + 1.14285 u^2 v^2 +
0.000884953 u v^3 + 1.39221*10^-6 v^4
The coefficients are not all minuscule, but all of the terms are quartic
which helps the objective function to stay close to zero over a sizeable
2-dimensional region around the solution point.
Check what the objective function looks like in this region by generating a
contour plot.
ContourPlot[objsoln[u, v], {u, -0.1, 0.1}, {v, -0.1, 0.1}]
which shows that the (u,0) and (0,v) directions are the "flattest".
Anyway, those are the sorts of techniques that I would use to get a "toe
hold" on the solution space around a particular found solution point.
--
Stephen Luttrell
West Malvern, UK
Prev by Date:
**Re: Can I solve this system of nonlinear equations?**
Next by Date:
**help with speed up**
Previous by thread:
**Re: Can I solve this system of nonlinear equations?**
Next by thread:
**Re: Can I solve this system of nonlinear equations?**
| |