newton

%Newton
clear all
close all
format long g

eps=1e-12;
nmax=150;

VT=26e-3;
R=1e3;
Is=1e-15;
E=5;
w=2*pi*1e3;
ns=2*4;
time=(2*pi/w)/ns;
vs=sqrt(2)*E*sin(w*time);
xstart=sqrt(2)*E/2;
%xstart=6;
x(1)=xstart;%initial guess
%f_cur=f(x(1));
%df_cur=df_dx(x(1));
f_cur=x(1)-R*Is*(exp((vs-x(1))/VT)-1); %initial fct
df_cur=1+(R*Is/VT)*(exp((vs-x(1))/VT)); %initial fct derivative

ncount=1;
while abs(f_cur)>eps & ncount<nmax
    ncount=ncount+1;
    x(ncount)=x(ncount-1)-inv(df_cur)*f_cur;
    %f_cur=f(x(ncount));
    %df_cur=df_dx(x(ncount);
    f_cur=x(ncount)-R*Is*(exp((vs-x(ncount))/VT)-1);
    df_cur=1+(R*Is/VT)*(exp((vs-x(ncount))/VT));
    hst(ncount-1,1)=x(ncount);
    hst(ncount-1,2)=f_cur;
end

figure
semilogy(hst(:,1),abs(hst(:,2)),'b*')
xlabel("Root")
ylabel("Residual")
title("Newton Raphson method")
grid on


disp("Approximated root")
disp(x(ncount))
disp("Final residual on f")
disp(abs(f_cur))
%disp("Final error on x")
%disp(abs(a(ncount)-b(ncount)))
disp("Number of iterations")
disp(ncount)
disp("Time instant [s]")
disp(time)
disp("Voltage [V]")
disp(x(ncount))