The Newton Raphson Method

% Author - Ambarish Prashant Chandurkar
%The Newton Raphson Method
clc;
close all;
clear all;
syms x;
f=x*exp(x)-1; %Enter the Function here
g=diff(f); %The Derivative of the Function
n=input('Enter the number of decimal places:');
epsilon = 5*10^-(n+1)
x0 = input('Enter the intial approximation:');
for i=1:100
     f0=vpa(subs(f,x,x0)); %Calculating the value of function at x0
     f0_der=vpa(subs(g,x,x0)); %Calculating the value of function derivative at x0
  y=x0-f0/f0_der; % The Formula
err=abs(y-x0);
if err<epsilon %checking the amount of error at each iteration
break
end
x0=y;
end
y = y - rem(y,10^-n); %Displaying upto required decimal places
fprintf('The Root is : %f \n',y);
fprintf('No. of Iterations : %d\n',i);

[Torsten]

Why did you change your former question ? Don't you think this is impolite against William Rose who tried to answer it ?

Answer 1

[William Rose]

@Bruce,

It appears that your x is a 2-coulmn array where column1 = theta, column 2=phi. I assume you have a function pendulum(), which returns a 1-by-2 vector whose values are .

Does the code you supplied run without error? The code you supplied appears to implement a forward Euler method. You will need to change what happens inside the loop for i=1:n,..., end, if you want to implement the midpoint method.

for i=1:n
    thalf = t(i) + dt/2; t(i+1)=t(i)+dt;
    xhalf = x(i,:) + (dt/2)*pendulum(thalf,xhalf);
    %The equaiton above is an implicit equaiton, because xhalf appears on
    %both sides.  You will have to do some analysis of the derivative
    %function inside pendulum() to solve this.  
    x(i+1,:) = 2*xhalf-x(i,:);
end

That gets you started, and you can finish the rest.