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: > > > > > > > > ((a-L)^2/a)*(2/(a+b))-p=0 > > > > ((b-U)^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, > > > > ((a-L)^2/a)*(2/(a+b))-p=0 > > ((b-U)^2/b)*(2/(a+b))-p=0 > > > > I get residuals > > > > 2.34718E-12 > > 1.9082E-16 > > > > 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 20-digit 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~ > 96738625-0.00067586120557976235658*#i,1.0006754048796738625-0.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]] > > Weierstrass-Durand-Kerner 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.
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.num-analysis. Thanks!