
Re: Terrible algebraic system
Posted:
Mar 7, 2012 3:36 PM


deltaquattro@gmail.com schrieb: > > Il giorno venerdì 24 febbraio 2012 17:21:59 UTC+1, Waldek Hebisch ha > scritto: > > > Hmmm! Tried following this approach, i.e. solving f(b)=0 with > > > Netwon's and then plugging b in a=g(b), but it doesn't work. The > > > solution doesn't satisfy the initial system. Try for example > > > > > > L=1 > > > U=1 > > > p=0.0013498980316301 > > > > > > I get > > > > > > a=1 > > > b=1 > > > > > > which clearly is wrong, i.e., a and b don't solve the initial > > > system, even though f(b)=f(1)=0 and a=g(b)=g(1)=1. I rewrite the > > > initial system here: > > > > > > ((aL)^2/a)*(2/(a+b))p=0 > > > ((bU)^2/b)*(2/(a+b))p=0 > > > > > > > In first message you wrote: > > > > > with a>L>0, b>U>0. > > > > Now you put L < 0, which violates this condition. > > > > > Am I doing something wrong, or is FriCAS/Derive solution wrong? > > > > FriCAS solution is "generic" in the sense that: > > > > 1) All solutions of your system are contained in FriCAS > > solution > > 2) It solves your system if there is no division by zero when > > plugging FriCAS solution into your system > > 3) Unless there are some special relation between L and U > > there will be no division by zero when plugging FriCAS > > solution into your system. > > > > In case of your system special relations are: L = 0, U = 0, > > L + U = 0. If no of the above relations hold then you > > need not worry about division by zero (in particular > > no worry when U>0 and L>0). When one of relations holds, > > you need to decide if you want to you your original > > system (that is throw out (a, b) which cause division > > by zero), or you want to simplify it and solve new, > > simpler system. > > > > BTW: Since a = 1, b = 1 is the only real solution that > > FriCAS found, this means that your orignal system even > > after simplification has no real solutions. > > > > Ok, tried again: I am using the formulation where a>L>0, b>U>0. Please > forget about the other! > > L = 1 > U = 1 > p = 0.022750132 > > With Newton's method, I solve > > 2 * b ^ 4 * (p  1) * (p  2) + 2 * b ^ 3 * (L * p * (p  2)  U * (p > ^ 2  8 * p + 8)) + b ^ 2 * (L ^ 2 * p ^ 2 + 8 * L * U * p + U ^ 2 * > (p ^ 2  16 * p + 24))  4 * b * U ^ 2 * (L * p + 2 * U * (2  p)) + 2 > * U ^ 4 * (2  p) = 0 > > for b, starting from b=U. I get b = 0.859373503 (!), which is not > bigger than U. I plug this result in > > a = (2 * b ^ 3 * (p  1) * (p  2) + 2 * b ^ 2 * (L * p * (p  2)  U > * (p ^ 2  8 * p + 8)) + b * (L ^ 2 * p ^ 2 + 8 * L * U * p + 4 * U ^ > 2 * (5  3 * p))  4 * U ^ 2 * (L * p + U * (2  p))) / (U ^ 2 * p * > (p  2)) > > getting a = 1.163638391, which is not smaller than L. Trying to plug > these solutions back in the original system, > > ((aL)^2/a)*(2/(a+b))p=0 > ((bU)^2/b)*(2/(a+b))p=0 > > I get residuals > > 2.34718E12 > 1.9082E16 > > which show that the system is indeed solved with excellent accuracy! > However, my constraints > > a>L > b>U > > are violated. Any idea how to "enforce" these constraints in the > problem solution? My guess: I should find all 4 solutions of the > equation for b by Newton's method, plug all of them one at a time in > the equation and discard all but the ones which satisfy constraints. > How can I find all the solutions with Newton's method? I usually only > use it to find a single zero. Can you help me? Thanks, >
Here are 20digit solutions for the cases you mentioned:
[L:=1,U:=1,p:=0.0013498980316301]
[[1,1],[1.0006754048796738625+0.00067586120557976235658*#i,1.0~ 006754048796738625+0.00067586120557976235658*#i],[1.00067540487~ 967386250.00067586120557976235658*#i,1.00067540487967386250.00~ 067586120557976235658*#i]]
[L:=1,U:=1,p:=0.022750132]
[[1.1776225370286649805,1.1776225370286649805],[0.85937350255160~ 052169,1.1636383912592832729],[1.1636383912592832729,0.859373502~ 55160052169],[0.86893696171357480864,0.86893696171357480864]]
WeierstrassDurandKerner iteration is related to Newton's method and produces all zeros of a polynomial in parallel. Complex arithmetic is needed for this since roots may be complex even for real polynomials.
Martin.

