
Re: Terrible algebraic system
Posted:
Mar 8, 2012 4:26 AM


Il giorno mercoledì 7 marzo 2012 21:36:02 UTC+1, clicl...@freenet.de ha scritto: > > > > 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.
Hi Martin!
Excellent :)your result confirms my idea that probably the best way to deal with this beast, is to compute all roots and then find the only one which satisfies a>L>0, b>U>0, i.e., [1.1776225370286649805,1.1776225370286649805]. However, I need to implement this in a code, since I must solve the system automatically for arbitrary inputs. Unluckily, the rest of the code is in VBA :( so I have no access to complex arithmetics. I could try to implement a class for that...or rewrite all the rest of the code in some different language. Also, I'm not familiar with this WDK method. Anyway, I guess this is now mainly a numerical mathematics topic, so I will also ask for suggestions on sci.math.numanalysis. Thanks!
Sergio

