% Pendulum BVP % % Dirichlet boundary conditions % global T alpha beta % T=2*pi; alpha=0.7; beta=0.7; % %hold off mm=2^3; hh=[]; tol=1e-8; %%%%%%%%%%%%%%%%%%% % computation of a reference solution % we will use it in place of the exact solution % we use the difference method with a very % fine grid, i.e. h=1/(2^15), to compute the reference solution % we use a smaller tolerance, tol/1000, in the Newton iteration, in the % computation of the reference solution %%%%%%%%%%%% changeinval=0; display('choose the initial guess for the Newton iteration') display('different initial guesses lead to different solutions of G(theta)=0') display('"1" is 0.7+sin(x/2), "2" is 0.7cos(x), ">2" is 0.7cos(x)+0.5sin(x)') changeinval=input('1,2,3,?') mm=mm*2^12; h=T/(mm); m=mm-1; xx=[h:h:T-h]'; theta=zeros(m,1); if (changeinval == 1 ) theta0=0.7*ones(m,1)+sin(xx/2);% else if (changeinval == 2) theta0=0.7*cos(xx); else theta0=0.7*cos(xx)+0.5*sin(xx); end end % theta=theta0; maxNit=100; % % Newton iteration to solve G(theta)=0 % for k=1:maxNit thetaold=theta; J=JacG(theta); g=G(theta); theta=theta-J\g; % difference between two successive iterations as stopping criterion nr=norm(theta-thetaold)/norm(theta0) if (nr