0

分解式を使用して、二次方程式を解くためのコードをMATLABに実装しています。

ここに画像の説明を入力してください

コードは次のとおりです。

clear all 
format short 
a=1; b=30000000.001; c=1/4; 
rdelta=sqrt(b^2-4*a*c); 
x1=(-b+rdelta)/(2*a);
x2=(-b-rdelta)/(2*a);
fprintf(' Roots of the polynomial %5.3f x^2 + %5.3f x+%5.3f \n',a,b,c)  
fprintf ('x1= %e\n',x1)
fprintf ('x2= %e\n\n',x2)
valor_real_x1= -8.3333e-009;
valor_real_x2= -2.6844e+007;

error_abs_x1 = abs (valor_real_x1-x1);
error_abs_x2 = abs (valor_real_x2-x2);

error_rel_x1 = abs (error_abs_x1/valor_real_x1);
error_rel_x2 = abs (error_abs_x2/valor_real_x2);

fprintf(' absolute_errorx1 = |real value - obtained value| = |%e - %e| = %e \n',valor_real_x1,x1,error_abs_x1)  
fprintf(' absolute_errorx2 = |real value - obtained value| = |%e - %e| = %e \n\n',valor_real_x2,x2,error_abs_x2) 

fprintf(' relative error_x1 = |absolut error / real value| = |%e / %e| = %e \n',error_abs_x1,valor_real_x1,error_rel_x1 ) 
 fprintf(' relative_error_x2 = |absolut error / real value|  = |%e / %e| = %e \n',error_abs_x2,valor_real_x2,error_rel_x2)

私が抱えている問題は、正確な解が得られることです。つまり、値a = 1、b = 30000000,001 c = 1/4の場合、根の値は次のようになります。

Roots of the polynomial 1.000 x^2 + 30000000.001 x+0.250 
 x1= -9.313226e-009
 x2= -3.000000e+007

多項式の根の正確な値は次のとおりです。

x1= -8.3333e-009
x2= -2.6844e+007

これにより、計算の絶対精度と相対精度に次のエラーが発生します。

 absolute_errorx1 = |real value - obtained value| = |-8.333300e-009 - -9.313226e-009| = 9.799257e-010 
 absolute_errorx2 = |real value - obtained value| = |-2.684400e+007 - -3.000000e+007| = 3.156000e+006 

 relative error_x1 = |absolut error / real value| = |9.799257e-010 / -8.333300e-009| = 1.175916e-001 
 relative_error_x2 = |absolut error / real value|  = |3.156000e+006 / -2.684400e+007| = 1.175682e-001

私の質問は次のとおりです。二次方程式の根を取得するための最適な方法はありますか？つまり、期待される解と結果の解の間の相対誤差を減らすためにコードを変更できますか？

4

4 に答える 4

7

この場合、二次方程式を直接使用すると、非常に類似した大きさの2つの値を減算すると、数値の精度が大幅に低下します。これは、式が

sqrt(b*b - 4*a*c)

bとほぼ同じです。したがって、これら2つの根の1つだけを使用する必要があります。一方は、2つの非常に近い値を減算する必要がなく、もう一方の根には、（たとえば）2次の根の積がc/aであるという事実を使用できます。ギャップを埋めさせてあげましょう。

于 2012-04-04T13:53:26.947 に答える
3

なぜこれは数値解析のファーストクラスからの宿題の問題のように聞こえますか？

幼い頃から久しぶりですが、思い出すとコツがあります。とにかく、あなたは間違っています。その多項式の真のルーツは

solve('x^2 + 30000000.001*x + 0.25')
ans =
          -30000000.000999991666666666944442
 -0.0000000083333333330555578703796293981491

根はここでどれくらいうまくいきますか？

p = [1 30000000.001 1/4];
format long g
roots(p)
ans =
             -30000000.001
     -8.33333333305556e-09

それは実際にはかなり良いようです。HPFはどのように機能しますか？

DefaultNumberOfDigits 64
a = hpf(1);
b = hpf('30000000.001');
c = hpf('0.25');

r1 = (-b + sqrt(b*b - 4*a*c))/(2*a)
r1 =
-0.000000008333333333055557870379629398149125529835186899898569329967

r2 = (-b - sqrt(b*b - 4*a*c))/(2*a)
r2 =
-30000000.000999991666666666944442129620370601850874470164813100101

はい、HPFも十分に機能します。

では、倍精度の数値と標準の数式を使用するとどうなりますか？ええ、クラポラが到着します。

a = 1;
b = 30000000.001;
c = 0.25;

>> r1 = (-b + sqrt(b*b - 4*a*c))/(2*a)
r1 =
     -7.45058059692383e-09

>> r2 = (-b - sqrt(b*b - 4*a*c))/(2*a)
r2 =
             -30000000.001

繰り返しになりますが、大規模な減算キャンセルは結果を食いつぶします。（それがあなたが最後の質問で抱えていた問題だったことを思い出しているようです。）

使用できるトリックがあります。ゼロに近い解ではなく、大規模な解が適切に推定されていることを確認してください。では、2次方程式を使用してfliplr（p）の根を解いた場合はどうなりますか？これはどのようにあなたの問題を解決しますか？あなたがそれをするとき、どんな変換が暗黙のうちに行われますか？（申し訳ありませんが、宿題はしません。とにかく、上記で十分なヒントになったと思います。）

于 2012-04-04T13:46:33.573 に答える
1

私はあなたの「本当の」値が間違っているかもしれないと思います（または多分それは正確なものです...私は知りません）

a*(valor_real_x1^2)+b*(valor_real_x1)+c

ans =

   9.9999e-07

a*(valor_real_x2^2)+b*(valor_real_x2)+c

ans =

  -8.4720e+13
于 2012-04-04T13:41:03.260 に答える
1

この問題の良い式：

var q = sqrt(c*a)/b;
var f = .5 + .5 *sqrt(1-4*q*q);
var x1=-b*f/a;
var x2=-c/(f*b);
于 2015-11-03T15:44:03.690 に答える