mit sehr kleinen Koeffizienten in Matlab eine quadratische Gleichung lösen

Question

ich einen Code in Matlab Implementierung quadratische Gleichungen zu lösen, die resolvent Formel:

Heres der Code:

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) 

Das Problem, das ich habe ist, dass es mir eine exakte Lösung, dh für Werte a = 1, b = 30000000.001 c = 1/4, die Werte der Wurzeln gibt:

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

Zu wissen, dass der genaue Wert der Wurzeln des Polynoms ist:

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

die mir die folgenden Fehler in der absoluten und relativen Genauigkeit der Berechnungen ergibt:

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 

Meine Frage ist, : Gibt es eine optimale Methode, um die Wurzeln einer quadratischen Gleichung zu erhalten? Dh ich kann Änderungen an meinem Code vornehmen, um den relativen Fehler zwischen der erwarteten Lösung und der resultierenden Lösung zu reduzieren?

Answer 1

Die Verwendung der quadratischen Formel direkt in diesen Fällen führt zu einem großen Verlust der numerischen Genauigkeit durch Subtrahieren von zwei Werten von sehr ähnliche Größenordnung.Dies liegt daran, dass der Ausdruck

fast identisch ist mit b. Sie sollten also nur eine dieser beiden Wurzeln verwenden, die keine Subtraktion von zwei sehr nahen Werten erfordert, und für die andere Wurzel können Sie zum Beispiel die Tatsache verwenden, dass das Produkt von Wurzeln einer Quadra- tik c/a ist. Ich lasse dich die Lücken füllen.

Answer 2

Ich denke, Ihre „echten“ Werte falsch sein könnten (oder vielleicht ist es eine Präzisions Sache ... ich weiß nicht)

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 

Answer 3

Warum klingt das wie ein Hausaufgabe Problem von einer ersten Klasse in der numerischen Analyse?

Es ist schon eine Weile her, seit ich so jung war, aber soweit ich mich erinnere gibt es einen Trick. Wie auch immer, du liegst falsch. Die wahren Wurzeln dieses Polynoms sind

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

Wie gut tut Wurzeln hier?

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

Das scheint eigentlich ziemlich gut. Wie macht 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 

Ja, HPF funktioniert auch gut genug.

Was passiert also, wenn Sie doppelt genaue Zahlen und die Standardformel verwenden? Ja, Crapola kommt an.

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 

Noch einmal, massive subtraktive Aufhebung frißt das Ergebnis auf. (Ich erinnere mich, dass das das Problem war, das Sie in Ihrer letzten Frage hatten.)

Es gibt einen Trick, den Sie verwenden können. Sehen Sie, dass die große Lösung gut geschätzt wurde, nur nicht die nahe Null. Was passiert also, wenn Sie mit der quadratischen Formel nach den Wurzeln von fliplr (p) suchen? Wie löst das dein Problem? Welche Transformation wird implizit vorgenommen, wenn Sie das tun? (Tut mir leid, aber ich werde deine Hausaufgaben nicht machen. Ich denke, das obige war sowieso genug eines Hinweises.)

Answer 4

Eine schöne Formel für dieses Problem:

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);