%
% equazione della catenaria risolta con il metodo di Newton
%
% soluzione di riferimento
clear all
close all 
a = 2;
mrif = 1001;
m = mrif;
x = linspace(-1,1,m)';
h = 2/(m-1);
A = toeplitz(sparse([1,1],[1,2],[-2/h^2,1/h^2],1,m));
A(1,1:2) = [1/h^2,0];
A(m,m-1:m) = [0,1/h^2];
B = toeplitz(sparse(1,2,-1/(2*h),1,m),sparse(1,2,1/(2*h),1,m));
B(1,2) = 0;
B(m,m-1) = 0;
F = @(u) A*u-a*sqrt(1+(B*u).^2)-[1/h^2-a;zeros(m-2,1);1/h^2-a];
J = @(u) A-a*spdiags(B*u./sqrt(1+(B*u).^2),0,m,m)*B;
tol = h^2/100;
u0 = ones(m,1);
u = u0;
errest = -J(u)\F(u);
while (norm(errest,inf) > tol)
  u = u+errest;
  errest = -J(u)\F(u);
end
urif = u;
mrange = [10,20,50,100];
counter = 0;
for m = mrange+1
  counter = counter+1;
  x = linspace(-1,1,m)';
  h = 2/(m-1);
  A = toeplitz(sparse([1,1],[1,2],[-2/h^2,1/h^2],1,m));
  A(1,1:2) = [1/h^2,0];
  A(m,m-1:m) = [0,1/h^2];
  B = toeplitz(sparse(1,2,-1/(2*h),1,m),sparse(1,2,1/(2*h),1,m));
  B(1,2) = 0;
  B(m,m-1) = 0;
  F = @(u) A*u-a*sqrt(1+(B*u).^2)-[1/h^2-a;zeros(m-2,1);1/h^2-a];
  J = @(u) A-a*spdiags(B*u./sqrt(1+(B*u).^2),0,m,m)*B;
  tol = h^2/100;
  u0 = ones(m,1);
  u = u0;
  errest = -J(u)\F(u);
  while (norm(errest,inf) > tol)
    u = u+errest;
    errest = -J(u)\F(u);
  end
  errorinf(counter) = norm(u-urif(1:(mrif-1)/(m-1):mrif),Inf);
end
loglog(mrange,errorinf,'*',mrange,errorinf(end)*(mrange/mrange(end)).^(-2),'k')

