15
2020
07

解一元二次方程

(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方以内,需要迭代几次,可以试试。



« 上一篇下一篇 »

发表评论:

◎欢迎参与讨论,请在这里发表您的看法、交流您的观点。