マラー法による3次方程式の解の数値解法

カルダノの方法によって、解は
=-2.637834253
=1.6810829-3.0504302i
=-1.6810829‐ +3.0504302i

マラー法では割線法に比べて解が高速に収束する。

clear; close all;

buffer=20;
x=zeros(1,buffer);
y=zeros(1,buffer);
h=zeros(1,buffer);
g=zeros(1,buffer);
l=zeros(1,buffer);

f = @(x) x^3+6*x^2+21*x+32;
f2 = @(x) x.^3+6*x.^2+21*x+32;

x(1)=0;x(2)=-1;x(3)=-2;
y(1)=f(x(1));y(2)=f(x(2));y(3)=f(x(3));

figure();
t=-4:.1:0;
plot(t,f2(t))
hold on;
plot(t,zeros(size(t)));
plot(x(1),0,'ro');
plot(x(2),0,'ro');
plot(x(3),0,'ro');
text(x(1),-2,'x(1)');
text(x(2),-2,'x(2)');
text(x(3),-2,'x(3)');

res=1e4;
delta=1e-8;

k=3;
while res>delta
    h(k)=x(k)-x(k-1);
    l(k)=h(k)/(x(k-1)-x(k-2));
    g(k)=(1+2*l(k))*(y(k)-y(k-1)) ...
        -(l(k))^2*(y(k-1)-y(k-2));
    
    D=(g(k))^2-4*y(k)*(1+l(k))*l(k)...
        *(y(k)-y(k-1)-l(k)*(y(k-1)-y(k-2)));
    denominator1 = g(k)+real(sqrt(D));
    denominator2 = g(k)-real(sqrt(D));
    if abs(denominator1)>=abs(denominator2) 
       denominator = denominator1;
    else
       denominator = denominator2;        
    end
    l(k+1)=-2*y(k)*(1+l(k))/denominator;
    
    x(k+1)=x(k)+h(k)*l(k+1);
    plot(x(k+1),0,'bx');
    txt = sprintf('x(%d)',k+1);
    text(x(k+1),5*randn(),txt);

    y(k+1)=f(x(k+1));
    h(k+1)=h(k)*l(k+1);    
   
    res = abs(f(x(k+1)));
    disp(res);
    k = k+1;
end

fprintf('ans =\n %d\n',x(k));

計算結果

ans =
 -2.637834e+00