why isn't my newton rhapson method giving me a quadratic conversion rate?

The code is running just fine, but the value konv(end)/(konv(end-1)^2) intended to show (e(i+1)/(e(i)^2)) gives an answer of 9.9952e+07, which seems to indicate that the conversion rate isn't quadratic. Is there something obviously wrong with the code? Help appreciated.

%constants

d=0.007;
dI=1.219;
Rfi=1.76*(10^(-4));
Rf0=Rfi;
hs=356;
ht=hs;
kw=60;
dTm=29.6;
Q=801368;
Aex=64.15;

a=((1/ht)+Rfi)/dI;
b=log(d/dI)/(2*kw);
c=d*(a+b)+Rf0+(1/hs);

A=((1/ht)+Rfi)/dI;
B=(log(d/dI)+1)/2*kw;

f=c*(Q/dTm)-Aex;

df=(Q/dTm)*(A+B);

%d1-do ska va mindre än tol
%i slutet av iterationen
%så man kan sätta while

d=0.007;
d0=d;
d1=10;
abs(d1-d0);

iteration=0;

konv=[0]

while abs(d1-d0)>10^(-8);
    
   a=((1/ht)+Rfi)/dI;
    b=log(d/dI)/(2*kw);
    c=d*(a+b)+Rf0+(1/hs);

    f=c*(Q/dTm)-Aex;

    A=((1/ht)+Rfi)/dI;
    B=(log(d/dI)+1)/2*kw;

    df=(Q/dTm)*(A+B);
    d0=d;
    d1=d-f/df;
    d=d1;
    iteration=iteration+1;   
    konv=[konv, abs(d1-d0)];
    
    
  
end

konv(end)/(konv(end-1)^2)

iteration
d

Answer 1

[Torsten]

If

b=log(d/dI)/(2*kw);

instead of

b=log(d/dI)/2*kw;

as was first written in your previous code, then B must also be

B=(log(d/dI)+1)/(2*kw);

instead of

B=(log(d/dI)+1)/2*kw;

In other words: df is wrong because B is wrong.

[Linnea]

it gives me the same answer for d though

[Torsten]

Yes, but did you see the number of iterations needed ? To reach convergence, the derivative needn't be exact, but as a consequence, convergence can become extremely slow.

[Linnea]

it was around 20 000, how do i make it converge faster?

[Torsten]

it was around 20 000, how do i make it converge faster?

By supplying a correct derivative df as I did. The number of iterations until convergence are 4 if you do so.

[Linnea]

could you please repeat what exactly you did? i changed B to B=(log(d/dI)+1)/(2*kw); but it's still giving me 20 000 iterations.

[Torsten]

Replace

%B=(log(d/dI)+1)/2*kw;

by

%B=(log(d/dI)+1)/(2*kw);

in the above code you posted. Nothing else is necessary.

%constants
d=0.007;
dI=1.219;
Rfi=1.76*(10^(-4));
Rf0=Rfi;
hs=356;
ht=hs;
kw=60;
dTm=29.6;
Q=801368;
Aex=64.15;
a=((1/ht)+Rfi)/dI;
b=log(d/dI)/(2*kw);
c=d*(a+b)+Rf0+(1/hs);
A=((1/ht)+Rfi)/dI;
B=(log(d/dI)+1)/2*kw;
f=c*(Q/dTm)-Aex;
df=(Q/dTm)*(A+B);
%d1-do ska va mindre än tol
%i slutet av iterationen
%så man kan sätta while
d=0.007;
d0=d;
d1=10;
abs(d1-d0);
iteration=0;
konv=[];
while abs(d1-d0)>10^(-10);
    
   a=((1/ht)+Rfi)/dI;
    b=log(d/dI)/(2*kw);
    c=d*(a+b)+Rf0+(1/hs);
    f=c*(Q/dTm)-Aex;
    A=((1/ht)+Rfi)/dI;
    B=(log(d/dI)+1)/(2*kw);
    df=(Q/dTm)*(A+B);
    d0=d;
    d1=d-f/df;
    d=d1;
    iteration=iteration+1;   
    konv=[konv, abs(d1-d0)];
    
    
  
end
for i = 2:numel(konv)
  konv(i)/konv(i-1)^2
end

ans = 17.0397

ans = 9.7667

ans = 9.1463

ans = 9.1980

iteration

iteration = 5

d

d = 0.0191