Reputation: 11
Newton’s method can be used to calculate the p-th root of a positive real number. Write a Matlab function myrootp.m which calculates a^(1/p) for any a> 0 and a positive integer p. The function should stop when abs((Xk^(p)/a)-1)<10^(-15).
I used the following code but I cannot find the error:
function root1 = root1( a,X0,p)
Xk=X0;
fd= p*Xk^(p-1);
f=(Xk^p)-a;
if a<0
then disp('Error:a must be positive and real')
else
for k=1:p
if abs(Xk^p/a-1)<10^(-15)
then return disp('Error');
else
Xk1=Xk-f/fd;
Xk=Xk1;
fd= p*Xk^(p-1);
f=(Xk^p)-a;
end
end
end
root1 = Xk;
return
end
Upvotes: 0
Views: 355
Reputation: 104503
I'd like to point out for those who are unfamiliar with this algorithm that this is finding the nth root of a number using Newton's method.
Back to your problem, there are two errors:
then
in MATLAB. Get rid of it.p
, and it's highly possible that the algorithm can converge before then, or the worst case require more iterations than p
to converge.Making these two changes:
function root1 = root1(a,X0,p)
Xk=X0;
fd= p*Xk^(p-1);
f=(Xk^p)-a;
if a<0
disp('Error:a must be positive and real')
else
while true
if abs(Xk^p/a-1)<10^(-15)
break;
else
Xk1=Xk-f/fd;
Xk=Xk1;
fd= p*Xk^(p-1);
f=(Xk^p)-a;
end
end
end
root1 = Xk;
To test this out, using an initial guess of 5 or X0=5
, and we want to find the square root of 10 (a = 10, p = 2
), we get:
>> format long g;
>> root1(10, 5, 2)
ans =
3.16227766016838
Comparing this with MATLAB's sqrt
function, we get:
>> format long g;
>> sqrt(10)
ans =
3.16227766016838
Upvotes: 1