(A+x)*(B+x) = w,已知A、B、w,求x。
=>A*B + (A+B)*x + x*x = w
=>x*x + (A+B)*x + A*B - w = 0
=>a = 1, b = A + B, c = A*B - w
然后带入一元二次方程的求根公式,就能得到两个根了。
本来以为就这么简单,然后就被坑了。
求根公式里边有个delta=b^2-4ac
当b*b>>4ac的时候(>>表示远大于),计算机认为,b^2-4ac=b^2,导致一个根是0。
可以用泰勒级数求出近似解。
设d=delta
√d = √(b*b-4ac) = b * √(1 - 4ac/b^2)
设x = 4ac/b^2,转化为b*√(1-x),泰勒级数展开根号下的内容
b*(1 - (1/2)x - (1/4)/(2!)x^2 - (3/8)/(3!)x^3-......);取前两项,误差量级为x^3,取前三项,误差量级为x^4,x很小,所以取前两项就足够精确了。
b - b * (0.5x + 0.125x^2)
r0 = -b + b - b * (0.5x + 0.125x^2) = -b * (0.5x + 0.125x^2)
r1 = -2b - r0
代码如下(typescript):
/** * 一元二次方程(A+x)*(B+x) = w * (A+x)*(B+x) = w * A*B + (A+B)*x + x*x = w * x*x + (A+B)*x + A*B - w = 0 * a = 1, b = A + B, c = A*B - w */ private oneVariableQuadraticEquation(A: number, B: number, w: number): number[] { const a: number = 1; const b: number = A + B; const c: number = A * B - w; const delta: number = b * b - 4 * a * c; // (A-B)^2+4*w一定大于0 const sqrDelta: number = Math.sqrt(delta); // 当b*b>>4ac时,算出来delta=b*b,解出来第一个结果为0. let r0: number = -b + sqrDelta; let r1: number = -b - sqrDelta; if (r0 === 0 && A * B !== w) { // b*b>>4ac导致的 // 设d=delta // √d = √(b*b-4ac) = b * √(1 - 4ac/b^2) // 设x = 4ac/b^2,转化为b*√(1-x),泰勒级数展开根号下的内容 // b*(1 - (1/2)x - (1/4)/(2!)x^2 - (3/8)/(3!)x^3-......);取前两项,误差量级为x^3,取前三项,误差量级为x^4,x很小,所以取前两项就足够精确了。 // b - b * (0.5x + 0.125x^2) // r0 = -b + b - b * (0.5x + 0.125x^2) = -b * (0.5x + 0.125x^2) // r1 = -2b - r0 const x: number = 4 * a * c / (b * b); r0 = -b * (0.5 * x + 0.125 * x * x); r1 = -2 * b - r0; } const a2: number = 1 / (2 * a); return [r0 * a2, r1 * a2]; }
直接用牛顿迭代法,收敛也很快。牛顿迭代法,解的误差控制在1e-30方以内,需要迭代几次,可以试试。
发表评论:
◎欢迎参与讨论,请在这里发表您的看法、交流您的观点。