Below you will find some Matlab scripts I wrote many years ago (2001-2002), for a postgraduate course on Numerical Analysis. These scripts solve twenty or so trivial problems using (and demonstrating) various numerical methods such as: Monte Carlo integration, uniform random deviates, Runge-Kutta with Bisection, Runge Kutta with Finite-Differences Newton-Raphson, Wilkinson deflation, Broyden-Charles and others.

Many of these scripts, are obsolete or deprecated, but it’s very easy to adapt them. Most of them can also be executed in Octave, with minor alterations…

```  % This Script finds the machine's precision using plain arithmetic. % % <- Dimitris B. Kalamaras ->   clear all; clc;   disp( 'Machine''s precision'); tol=1.00; e=1.00+tol; while ( e > 1.0000000 ) tol=tol/2.00; e=1.00+tol; end disp(' '); disp('TOL of the machine is: '); disp(tol);             % % Script for the approximation via Continuous Fraction Expansion of 2^sqrt(2) % % It seems that the continued fraction expansion is smaller and then bigger % and then smaller again and then bigger again and so on.... % % <- Kalamaras Dimitris ->   clc; echo off; clear all;   disp( 'Approximation via Continuous Fraction Expansion of 2^sqrt(2)');   n=input ('Enter Expansion Depth (1..40): ');   format long; r=zeros(n,1); q=zeros(n,1); a=zeros(n,1);   a(1)=vpa(2^sqrt(2)); % Enter another number if you wish...to approximate   for k=1:n q(k)=fix ( a(k) ); %finds the integer part of a(k) r(k)= a(k)-q(k) ; %finds the remainder a(k+1)= 1/ r(k); %inverts the remainder, which is a(k+1) end   for k=1:n number = help_function (1,k,q,r); if (double(a(1))-double(number))> 0 sprintf('Number No %d is %30.15f',k, number) sprintf('and is SMALLER than 2^sqrt(2) = %30.15f', a(1)) elseif (double(a(1))-double(number))< 0 sprintf('Number No %d is %30.15f',k,number) sprintf('and is BIGGER than 2^sqrt(2) = %30.15f', a(1)) else sprintf('This machine reached its precision limits: 2^sqrt(2) = %30.15f', number ) end end format;         function y = help_function (i,n,q,r) % % % The "help_function" is used for the computation of the Continued Fraction % % % This "help_function" uses the idea of recursion. % A function that calls itself is said to be recursive. % % <- Dimitris Kalamaras ->   % clc;   if nargin==0 error('Usage: help_function(i,n,q,r), where n integer'); end   if i < n y=q(i)+1/help_function(i+1,n,q,r); else y=q(n); end               % % This script uses the Newton Raphson method in order to find the solutions % of the equation : % cos(x)-exp(-x^2)=0 % % The scripts demands an estimation x0 and returns the approximate root. % % % <-Kalamaras Dimitris->   clc; clear all; format long; %x=zeros(140,1);   disp(' NEWTON-RAPHSON method the equation : '); disp(' cos(x)-exp(-x^2)=0 ');   x(1)=input('Please enter an estimation of the root x0 : '); repeats=input('Enter amount of repeats for the Newton-Raphson Method : '); tol=input('Enter desired precision :');   for k=1:repeats f=[cos(x(k))-exp( -x(k)*x(k) )]; df=[-sin(x(k))+2*x(k)*exp( -x(k)*x(k) )];   x(k+1)=x(k)-f/df ;   if abs( x(k+1)-x(k) )< tol root=x(k+1); sprintf('The root is : %20.10f found with %d repeats', root, k+1) break end end   if abs( x(k+1)-x(k) )> tol sprintf('Could not converge to a solution with the given precision...!') end                 % This "nr" function uses Newton Raphson method for finding the roots of a % given function f % % <- Kalamaras Dimitris -> % clear all; clc; format long;   x=zeros(40,1);   disp('Newton Raphson method for the equation cos(x)-exp(-0.5x)');   x(1)=input('Enter guess value x(1) :');   %f=[1 0 0 -1] %df=polyder(f) %r=roots(f)   for k=1:40 f=cos(x(k))-exp(-0.5*x(k)); df=-sin(x(k))-0.5*exp(-0.5*x(k));   x(k+1)=x(k)+f/df; if (x(k+1)-x(k))< 0.001 break end end   root=x(k+1); disp('The solution of cosx-exp(-0.5x)=0 in the given interval is :'); disp(root);             % % This script uses the Newton Raphson method in order to find the solutions % of the equation : % cos(x)-exp(-x^2)=0 % % The scripts demands an estimation x0 and returns the approximate root. % % % <-Kalamaras Dimitris->   clc; clear all; format long; %x=zeros(140,1);   disp(' NEWTON-RAPHSON with FINITE DIFFERENCES method for the equation : '); disp(' cos(x)-exp(-x^2)=0 ');   x(1)=input('Please enter an estimation of the root x0 : '); repeats=input('Enter amount of repeats for the Newton-Raphson Method : '); tol=input('Enter desired precision :'); h=0.005; dummy=input('Starting with h = 0.005. PRESS '); for i=1:10 for k=1:repeats f=[cos(x(k))-exp( -x(k)*x(k) )];   fh= [cos(x(k)-h)-exp( -(x(k)-h)^2 )];   df=(f-fh)/h;   x(k+1)=x(k)-f/df ;   if abs( x(k+1)-x(k) )< tol root=x(k+1); sprintf('The root is : %20.10f found with %d repeats', root, k+1) break end end   if abs( x(k+1)-x(k) )> tol sprintf('Could not converge to a solution with the given precision...!') end   h=h*2; sprintf('Now step is h= %6.5f. PRESS ', h) dummy=input (' ' ); end           % This script is an implemention of the Bisection (Bolzano) % method which finds a root of a function, manifestated in f_bis.m, % in a given space (a,b) (explicitly). % % <- Kalamaras Dimitris -> % clear all; clc;   format long;   disp('BISECTION method for the equation cos(x)-exp(-x^2)=0'); a=input('Enter left interval bound a : '); b=input('Enter right interval bound b : '); repeats=input('Enter amount of repeats : '); tol=input('Enter desired precision : '); x=zeros(41,1); x0=a;   x(1)=x0+((b-a)/2^1)*sign(f_bis(a))*sign(f_bis(x0));   for k=1:repeats x(k+1)=x(k)+((b-a)/2^(k+1))*sign(f_bis(a))*sign(f_bis(x(k)));   if abs(x(k+1)-x(k))< tol root=x(k+1); sprintf('The root in the given interval is : %20.15f in %d times', root, k+1) break ; end   end   if abs( x(k+1)-x(k) ) > tol sprintf('Could not converge to a solution with the given precision ') end       function y=f_bis(x) % function 'f_bisn' is the the function % used in the 'bisect' method % % <- Kalamaras Dimitris ->   global c;   if nargin==0 x=0:.1:10; end %if nargin==0 error('Usage: f_bis( x ) , where x whatever.'); end   y=cos(x)-exp(-x*x);   % with some matlab native noise!   %noise=-ran0*c; %noise=randn*c; %y=y+noise*y;             % % This script uses the Newton Raphson method in order to find the solutions % of the equation : % cos(x)-exp(-x^2)=0 % % The scripts demands an estimation x0 and returns the approximate root. % % % <-Kalamaras Dimitris-> clc;   x=zeros(140,1);   disp(' SECANT method for the equation :'); disp(' cos(x)-exp(-x^2)=0 ');   x(1)=input('Please enter a first value x0 : '); x(2)=input('Please enter a first value x1 : ');   repeats=input('Enter amount of repeats for the Secant Method : '); tol=input('Enter desired precision : ');   for k=2:repeats   fkp=[cos(x(k-1))-exp( -x(k-1)*x(k-1) )]; fk=[cos(x(k))-exp( -x(k)*x(k) )]; x(k+1)=x(k)- fk * ( x(k)-x(k-1))/(fk-fkp);   if abs(x(k+1)-x(k))< tol root=x(k+1); sprintf('The root is : %15.5f found with %d repeats', root, k+1) break; end end   if abs(x(k+1)-x(k))> tol break; sprintf('Could not converge to a solution with the given precision') end               % This script examines the behaviour of Newton Raphson and BISECTION methods % concerning convergence. It applies white GAUSSIAN noise with various values % of variance and mean value zero for the solution of the equation : % cos(x)-exp(-x^2)=0 % % <- Dimitris B. Kalamaras ->   clc; clear all; format long; disp(' Convergence of NEWTON-RAPHSON & BISECTION methods for the equation : '); disp(' cos(x)-exp(-x^2)=0 '); disp(' by applying white Gaussian Noise for various sigmas'); disp(' ');   disp('NEWTON RAPHSON: estimation x(1)=5'); x(1)=5; disp('NEWTON RAPHSON & BISECTION: repeats 300 max'); repeats=300; disp('NEWTON RAPHSON & BISECTION: precision 0.0001'); tol=0.0001; disp('BISECTION: left interval bound a = 2 '); a=2; disp('BISECTION: right interval bound b =5 '); b=5; s=0.00001; key=0; disp(' '); dummy=input('The first results will have no Gaussian Noise. PRESS ');   for i=1:10 % here starts the NEWTON RAPHSON algorithm... for k=1:repeats f=[cos(x(k))-exp( -x(k)*x(k) )]+key*gaussian(0,s); df=[-sin(x(k))+2*x(k)*exp( -x(k)*x(k) )];   x(k+1)=x(k)-f/df ;   if abs( x(k+1)-x(k) )< tol root=x(k+1); sprintf('NEWTON - RAPHSON : The root is %20.10f found with %d repeats', root, k+1) break end end   if abs( x(k+1)-x(k) )> tol sprintf('Could not converge to a solution with the given precision...!') end     % Here starts the BISECTION algorithm...   x0=a;   fa=[cos(a)-exp( -a*a )]+key*gaussian(0,s); fx0=[cos(x0)-exp( -x0*x0 )]+key*gaussian(0,s);   x(1)=x0+((b-a)/2^1)*sign(fa)*sign(fx0);   for k=1:repeats fa=[cos(a)-exp( -a*a )]+key*gaussian(0,s); fxk=[cos(x(k))-exp( -x(k)*x(k) )]+key*gaussian(0,s); x(k+1)=x(k)+( (b-a)/2^(k+1))*sign(fa)*sign(fxk) ;   if abs(x(k+1)-x(k))< tol root=x(k+1); sprintf('BISECTION : The root in the given interval is %20.10f in %d repeats', root, k+1) break ; end   end   s=s*3;   sprintf('Now Using Gaussian Noise with sigma s= %8.5f . PRESS ',s) dummy=input('');   clc; x(1)=5; key=1; end       disp('Apparently, the BISECTION method endures noise for all sigma''s');     function Y=gaussian ( mean, sigma) % % Function to return random variable distributed % according to Gaussian distribution with mean mean % and standard deviation sigma. % It uses the "ACCEPT-REJECT" method % % < Dimitris B. Kalamaras >   if nargin < 2 error('Usage : gaussian (mean, sigma) '); end;   ymin = mean - 4. * sigma; ymax = mean + 4. * sigma; Pymax = 1. / sqrt (2. * pi) / sigma; % this is pymax .... % Calculate random value uniformly distributed % in range ymin to ymax y = ymin + (ymax - ymin) * rand; %/ double (RANDMAX); % Calculate Py Py = exp (- (y - mean) * (y - mean) / 2. / sigma / sigma) / sqrt (2. * pi) / sigma; % Calculate random value uniformly distributed in range 0 to Pymax x = Pymax * rand ; %/ double (RANDMAX); % If x > Py reject value and recalculate if x > Py % disp(x), disp(Py), disp(y); Y=gaussian (mean, sigma); else Y=y; end             % % This script implements the WILKINSON deflation technique using Newton Raphson method to find % the solutions of the equation : % (x-1)^2*(x+2)*x=0 % % The scripts demands assumes x0=4 and returns the approximate roots over the domain of the function. % % % <-Kalamaras Dimitris-> %   clc; clear all; format long;   disp(' Implementation of the WILKINSON technique using NEWTON-RAPHSON method for the equation : '); disp(' (x-1)^2(x+2)x=0 '); disp('Estimation of the root: x0 = 4 ');   x(1)=4; repeats=input('Enter amount of repeats for the Newton-Raphson Method ( > 50 ): '); tol=input('Enter desired precision [ < 0.0000001 ]:');   %first loop - first root (x0=1 multiplicity 2) root_number=1; for k=1:repeats f=[(x(k)-1)^2 * (x(k)+2)*x(k) ]; df=[2*( x(k)-1 ) * (x(k)+2)* x(k) + (x(k)-1)^2 * x(k) + (x(k)-1)^2 * (x(k)+2)];   x(k+1)=x(k)-f/df ;   if abs( x(k+1)-x(k) )< tol root(root_number)=x(k+1); sprintf('The %1dst root is : %20.10f found with %d repeats',root_number, root(root_number), k+1) root_number=root_number+1; for i=1:k+1 % vector of repeats used in plot "against" command.... y(i)=i; end plot(y,x,'rx:'); clear x; x(1)=4; %second loop - second root (x0=1 multiplicity 2) for k=1:repeats f=[ (x(k)-1) * (x(k)+2)*x(k) ]; df=[ (x(k)+2)* x(k) + (x(k)-1) * x(k) + (x(k)-1) * (x(k)+2)];   x(k+1)=x(k)-f/df ;   if abs( x(k+1)-x(k) )< tol root(root_number)=x(k+1); sprintf('The %1dst root is : %20.10f found with %d repeats',root_number, root(root_number), k+1) root_number=root_number+1; clear y; for i=1:k+1 % vector of repeats used in plot "against" command.... y(i)=i; end hold on ; plot(y,x,'bx:'); clear x; x(1)=4; %third loop - third root (x0=0) for k=1:repeats f=[ (x(k)+2)*x(k) ]; df=[ x(k) + (x(k)+2)];   x(k+1)=x(k)-f/df ;   if abs( x(k+1)-x(k) )< tol root(root_number)=x(k+1); sprintf('The %1dst root is : %20.10f found with %d repeats',root_number, root(root_number), k+1) root_number=root_number+1; clear y; for i=1:k+1 % vector of repeats used in plot "against" command.... y(i)=i; end hold on;plot(y,x,'go:'); clear x; x(1)=4; %fourth loop - fourth root (x0=-2) for k=1:repeats f=[ (x(k)+2) ]; df=[1];   x(k+1)=x(k)-f/df ;   if abs( x(k+1)-x(k) )< tol root(root_number)=x(k+1); sprintf('The %1dst root is : %20.10f found with %d repeats',root_number, root(root_number), k+1) root_number=root_number+1; clear y; for i=1:k+1 % vector of repeats used in plot "against" command.... y(i)=i; end hold on;plot(y,x,'yd:'); text (10,3,'Red no deflation, Blue 1st,Green 2ond and Yellow 4th deflation.'); break end end   if abs( x(k+1)-x(k) )> tol sprintf('Could not converge to a solution with the given precision...!') end     break end end   if abs( x(k+1)-x(k) )> tol sprintf('Could not converge to a solution with the given precision...!') end     break end end   if abs( x(k+1)-x(k) )> tol sprintf('Could not converge to a solution with the given precision...!') end     break end end   if abs( x(k+1)-x(k) )> tol sprintf('Could not converge to a solution with the given precision...!') end             % % CONTOUR plots for the the system % f1= x.^2-y.^2-1 ; % f2= x.^2+y.^2-2; % % % <- Kalamaras Dimitris ->   function contr1(x1,x2);   if nargin==0 dx=1/8; [x y]=meshgrid(-3:dx:3); end;       disp( 'CONTOUR plots for the the system : '); disp( ' f1= x.^2-y.^2-1 '); disp( ' f2= x.^2+y.^2-2 ');       f1= x.^2-y.^2-1 ; f2= x.^2+y.^2-2; mesh(x,y,f1); hold on; mesh(x,y,f2); title('Mesh plot for f1=x^2-y^2-1 and f2=x^2+y^2-2');   figure surf(x,y,f1); hold on; surf(x,y,f2); title('Surf plot for f1=x^2-y^2-1 and f2=x^2+y^2-2');   figure contour(x,y,f1); title('Contour plot for f1=x^2-y^2-1');   figure contour(x,y,f2); title('Contour plot for f2=x^2+y^2-2');   figure contour(x,y,f1); hold on; contour(x,y,f2); title('Contour plot for f1=x^2-y^2-1 and f2=x^2+y^2-2');               % % This Matlab script approximates the solution of the system % {f1==x1^2 + x2^2 - 9 = 0 % {f2==x1 + x2 - 3 = 0 , where f1,f2 are R^2-->R, % using the Newton Method in 2 dimension: % Jf(xp)*Sp=-F(xp) % x(p+1)=x(p)+Sp % % % <-Kalamaras Dimitris-> %   clear all; clc;     disp('NEWTON 2d for the solution of the following system:'); disp(' f1(x1,x2)=x1^2 + x2^2 - 9'); disp(' f2(x1,x2)=x1 + x2 - 3'); disp(' with initial values (x1,x2)=(1,5)'); % f1(x1,x2)=x1^2-x2^2-9 % f2(x1,x2)=x1+x2-3     x1=1; %initial values x2=5;   max_repeats=input('Enter maximum amount of repeats for Newton 2d :');   %If you want another system, change the two next lines. fxr=[x1^2+x2^2-9 ; x1+x2-3]; %Functions' f1 & f2 values vector Jf=[2*x1 , 2*x2; 1 , 1 ]; % Jacobian of f1, f2 ref to x1,x2   for k=1:max_repeats   Sp=Jf\(-fxr); % Left division of Jfxr into -fxr uses LU factorization x1p=x1; x2p=x2; x1=x1+Sp(1); x2=x2+Sp(2); if abs(x1-x1p)<0.0001 & abs(x2-x2p)<0.0001 break; end;   fxr=[x1^2+x2^2-9 ; x1+x2-3]; %Functions' values NEW vector Jf=[2*x1 , 2*x2; 1 , 1 ]; % NEW Jacobian of f1, f2 ref to x,y end;   disp('Solution is:'); root=[x1;x2]; disp(root);   disp('iterations done : '); disp(k);   %toc             % This Matlab script approximates the solution of the system % {f1==x1^2 + x2^2 - 9 = 0 % {f2==x1 + x2 - 3 = 0 , where f1,f2 are R^2-->R, % % using the Broyden Charles Method which is a quasi-Newton in 2 dimensions: % % A(p)*Sp=-F(x(p)) % x(p+1)=x(p)+Sp % y(p)=F(x(p+1))-F(x(p)) % A(p+1)=A(p)+((y(p)-A(p)*Sp)*Sp')/Sp'*Sp % % Broyden-Charles approximates the Jacobian, thus avoiding the calculation, if it is possible, % of the derivatives. % % Use as initial matrix A, a non singular, ie A=I. A should be converging to Jacobian % at the end of the Broyden-Charles. % % <-Dimitris B. Kalamaras -> %     clear all; clc;   disp(' BROYDEN method for the following system:'); disp(' f1(x1,x2)=x1^2 + x2^2 - 9 '); disp(' f2(x1,x2)=x1 + x2 - 3 '); disp(' with initial values (x1,x2)=(1,5)'); % f1(x1,x2)=x1^2+x2^2-9 % f2(x1,x2)=x1+x2-3     x1=1; %initial values x2=5;   max_repeats=input('Enter maximum amount of repeats for Broyden :');   A=eye(2); %Change to whatever non singular matrix 2x2.     %If you want another system, change the next lines. fxp=[x1^2+x2^2-9 ; x1+x2-3]; %Functions' f1 & f2 values vector   for k=1:max_repeats   Sp=A\(-fxp); % Left division of Jfxr into -fxr uses LU factorization x1p=x1; x2p=x2; x1=x1+Sp(1); x2=x2+Sp(2); if abs(x1-x1p)<0.0001 & abs(x2-x2p)<0.0001 break; end;   fxp1=[x1^2-x2^2-9 ; x1+x2-3]; %Functions' values NEW vector   x2p=fxp1-fxp; A=A+((x2p-A*Sp)*Sp')/(Sp'*Sp);   fxp=fxp1; end;   disp('The Broyden Solution is:'); root=[x1; x2]; disp(root);   disp('iterations done : '); disp(k);   %toc           % This Matlab script approximates the solution of the system % {f1==x1^2 + x2^2 - 9 = 0 % {f2==x1 + x2 - 3 = 0 , where f1,f2 are R^2-->R, % using the quasi-Newton Method with finite differences in 2 dimensions: % Jf(xp)*Sp=-F(xp) where jacobian is approximated % x(p+1)=x(p)+Sp % % % <-Kalamaras Dimitris->   clear all; clc;   disp(' Quasi-NEWTON with finite differences'); disp(' for the solution of the following system:'); disp('f1(x1,x2)=x1^2-x2^2-9'); disp('f2(x1,x2)=x1+x2-3'); disp('with initial values (x1,x2)=(1,5)');     %if you want to solve another system , change it in the function 'f_qn' ....   x1=1; %initial values x2=5;   max_repeats=input('Enter maximum amount of repeats for Newton 2d f.diff. :'); h=input('Enter step h : ');   theta1_f1=(f_qn(1,x1+h,x2)-f_qn(1,x1,x2))/h; theta2_f1=(f_qn(1,x1,x2+h)-f_qn(1,x1,x2))/h; theta1_f2=(f_qn(2,x1+h,x2)-f_qn(2,x1,x2))/h; theta2_f2=(f_qn(2,x1,x2+h)-f_qn(2,x1,x2))/h;     fxr=[f_qn(1,x1,x2) ;f_qn(2,x1,x2) ]; %Functions' f1 & f2 values vector Jf=[theta1_f1 , theta2_f1; theta1_f2 , theta2_f2 ]; % Jacobian of f1, f2 ref to x1,x2   for k=1:max_repeats   Sp=Jf\(-fxr); % Left division of Jfxr into -fxr uses LU factorization x1p=x1; x2p=x2; x1=x1+Sp(1); x2=x2+Sp(2); if abs(x1-x1p)<0.0001& abs(x2-x2p)<0.0001 break; end;   theta1_f1=(f_qn(1,x1+h,x2)-f_qn(1,x1,x2))/h; theta2_f1=(f_qn(1,x1,x2+h)-f_qn(1,x1,x2))/h; theta1_f2=(f_qn(2,x1+h,x2)-f_qn(2,x1,x2))/h; theta2_f2=(f_qn(2,x1,x2+h)-f_qn(2,x1,x2))/h;   fxr=[f_qn(1,x1,x2) ;f_qn(2,x1,x2) ]; %Functions' values NEW vector Jf=[theta1_f1 , theta2_f1; theta1_f2 , theta2_f2 ]; % NEW Jacobian of f1, f2 ref to x1,x2   end;   disp('NEWTON finite differences found the Solution to be: '); root=[x1;x2]; disp(root); disp('iterations done : '); disp(k); %toc         % % This script calls the 'newton21' (Newton for systems) function to solve the equation % z^3-1-0 and reports if a given point 'goes to' which root of the equation. % % It constructs the corresponding fractal, which shows the areas % of convergence in the given space [-2,2], [-2i,2i]. % % % <- Dimitris B. Kalamaras ->     clc; clear all; tic;     disp( 'Generation of fractal from the solution of z^3-1=0 '); disp( ' using Newton for systems.');   i1=1; i2=1 ; i3=1; root1=[1;0]; root2=[-1/2;-sqrt(3)/2]; root3=[-1/2;sqrt(3)/2]; max_iterations=input('Enter number of repeats for Newton2D :'); for r=-2:0.025:2 %0.025 ==== 232 secs on celeron 333mhz for im=-2:0.025:2 k=newton21(r,im,max_iterations); disp('With initial value'); disp(r+im*i); disp('Newton2D converges to'); disp(k); if abs(k-root1) < [0.1;0.1] red(i1)=r+im*i; i1=i1+1; elseif abs(k-root2) < [0.1;0.1] green(i2)=r+im*i; i2=i2+1; elseif abs(k-root3)< [0.1;0.1] blue(i3)=r+im*i; i3=i3+1; end end end   figure; box off; hold on;   plot(blue,'b.'); plot(red,'r.'); plot(green,'g.'); title('Fractal, areas of convergence in [-2,2][-2i,2i] for z^3-1=0 via Newton in 2 dimensions'); xlabel('red=1, green=-1/2-sqrt(3)/2, blue=-1/2+sqrt(3)/2'); hold off; zoom on; toc;         function root=newton21 (a,b,max_iter) % Newton2 approximates the solution of the system {f1==0,f2==0}, where f1,f2 are R^2-->R, % using the Newton Method in 2 dimension: % Jf(xp)*Sp=-F(xp) % x(p+1)=x(p)+Sp   % % <-Kalamaras Dimitris->     %clear all; clc; if nargin<3 error('Usage newton21(a,b,max_iter) where a,b rational Re(z0) & Im(z0), max_iter integer'); end;   % for z^3-1=0 it is % z^3-1=0 % (x-y*i)^3-1=0 % x^3-3*i*x^2*y-3*x*y^2+i*y^3-1=0 % f1=Re(z^3-1)=x^3-3*x*y^2-1 % f2=Im(z^3-1)=y^3-3*x^2*y %tic     x=a; %initial values y=b;     %If you want another system, change the two next lines.   fxr=[x^3-3*x*y^2-1; y^3-3*x^2*y;]; %Functions' f1 & f2 values vector   Jf=[3*x^2-3*y^2, -6*x*y ; -6*x*y, 3*y^2-3*x^2]; % Jacobian of f1, f2 ref to x,y     for iteration=1:max_iter   Sp=Jf\(-fxr); % Left division of Jfxr into -fxr uses LU factorization xp=x; yp=y; x=x+Sp(1); y=y+Sp(2); if abs(x-xp)<0.0001 & abs(y-yp)<0.0001 break; end;   fxr=[x^3-3*x*y^2-1; y^3-3*x^2*y]; %Functions' values NEW vector   Jf=[3*x^2-3*y^2, -6*x*y ; -6*x*y, 3*y^2-3*x^2]; % NEW Jacobian of f1, f2 ref to x,y     end;   %disp('Solution is:'); Jf iteration root=[x;y];   %toc       % % This script calls the 'newton21_noise' (Newton for systems with gaussian noise) function to solve the equation % z^3-1-0 and reports if a given point 'goes to' which root of the equation. % % It constructs the corresponding fractal, which shows the areas % of convergence in the given space [-2,2], [-2i,2i]. % % % <- Dimitris B. Kalamaras -> %   clc; clear all; tic;     disp( 'Generation of fractal from the solution of z^3-1=0 '); disp( ' using Newton for systems with gaussian noise N(0,s).'); i1=1; i2=1 ; i3=1; root1=[1;0]; root2=[-1/2;-sqrt(3)/2]; root3=[-1/2;sqrt(3)/2]; max_iterations=input('Enter number of repeats for Newton2D :'); s=input('Enter value of parameter sigma for the gaussian noise : '); for r=-2:0.025:2 %0.025 ==== 232 secs on celeron 333mhz for im=-2:0.025:2 k=newton21_noise(r,im,max_iterations,s); %calls newton21 with noise to return a point... disp('With initial value'); disp(r+im*i); disp('Newton2D converges to'); disp(k); if abs(k-root1) < [0.1;0.1] red(i1)=r+im*i; i1=i1+1; elseif abs(k-root2) < [0.1;0.1] green(i2)=r+im*i; i2=i2+1; elseif abs(k-root3)< [0.1;0.1] blue(i3)=r+im*i; i3=i3+1; end end end   figure; box off; hold on;   plot(blue,'b.'); plot(red,'r.'); plot(green,'g.'); title('Fractal- areas of convergence in [-2,2][-2i,2i] for z^3-1=0 via Newton 2d with Noise s=0.1'); xlabel('red=1, green=-1/2-sqrt(3)/2, blue=-1/2+sqrt(3)/2'); hold off; zoom on; toc;             function root=newton21_noise (a,b,max_iter, s) % Newton2 approximates the solution of the system {f1==0,f2==0}, where f1,f2 are R^2-->R, % using the Newton Method in 2 dimension: % Jf(xp)*Sp=-F(xp) % x(p+1)=x(p)+Sp   % % <-Kalamaras Dimitris->     %clear all; clc; if nargin<3 error('Usage newton21_noise(a,b,max_iter,s) where a,b rational Re(z0) & Im(z0), max_iter integer, s: sigma'); end;   % for z^3-1=0 it is % z^3-1=0 % (x-y*i)^3-1=0 % x^3-3*i*x^2*y-3*x*y^2+i*y^3-1=0 % f1=Re(z^3-1)=x^3-3*x*y^2-1 % f2=Im(z^3-1)=y^3-3*x^2*y     x=a; %initial values y=b;     %If you want another system, change the two next lines.   fxr=[x^3-3*x*y^2-1+gaussian(0,s); y^3-3*x^2*y+gaussian(0,s);]; %Functions' f1 & f2 values vector   Jf=[3*x^2-3*y^2, -6*x*y ; -6*x*y, 3*y^2-3*x^2]; % Jacobian of f1, f2 ref to x,y     for iteration=1:max_iter   Sp=Jf\(-fxr); % Left division of Jfxr into -fxr uses LU factorization xp=x; yp=y; x=x+Sp(1); y=y+Sp(2); if abs(x-xp)<0.0001 & abs(y-yp)<0.0001 break; end;   fxr=[x^3-3*x*y^2-1; y^3-3*x^2*y]; %Functions' values NEW vector   Jf=[3*x^2-3*y^2, -6*x*y ; -6*x*y, 3*y^2-3*x^2]; % NEW Jacobian of f1, f2 ref to x,y   end;   %disp('Solution is:'); Jf iteration root=[x;y];           % % ITERATIVE METHODS % % JACOBI - GAUSS SEIDEL - S.O.R. METHODS FOR SOLVING LINEAR NXN SYSTEMS % % % <-Dimitris B. Kalamaras -> %     clc; clear all;     disp( 'LINEAR SYSTEMS SOLVE METHODS (iterative)'); disp(' JACOBI - GAUSS SEIDEL - S.O.R. methods for a NxN linear system '); disp('-------------------------------------------------------------------------------- ');   disp( ' System''s dimension, N '); n=4; disp(n); disp( 'Matrix A, '); A(1,1)=4; A(1,2)=-1; A(1,3)=-1; A(1,4)=0; A(2,1)=-1; A(2,2)=4; A(2,3)=0; A(2,4)=-1; A(3,1)=-1; A(3,2)=0; A(3,3)=4; A(3,4)=-1; A(4,1)=0; A(4,2)=-1; A(4,3)=-1; A(4,4)=4;   disp(A);   disp( 'Matrix B,'); B(1)=0; B(2)=0; B(3)=1000; B(4)=1000; disp(B); disp( 'Matrix X,'); X(1)=1; X(2)=1; X(3)=1 ; X(4)=1; disp(X);   max_repeats=15; disp('Maximum repeats for all three methods, '); disp(max_repeats); OMEGA=1.062; disp('relaxation parameter OMEGA for S.O.R. : '); disp(OMEGA); TOL=0.0001; disp('precision TOL : '); disp(TOL); dummie=input('Press to Continue>');   %main JACOBI-GAUSS SEIDEL- S.O.R. algorithms (repeated for algorithmical elegance)   for method=1:3   for repeats=1:max_repeats   for k=1:n previous_X(k)=X(k); end   for i=1:n   sum1=0; sum2=0;   if i~=1 for j=1:i-1 sum1=sum1+A(i,j)*X(j); end end;   if i~=n for j=i+1:n sum2=sum2+A(i,j)*X(j); end end;   if method==1 Z(i)=(1/A(i,i))*(B(i)-sum1-sum2); elseif method==2 X(i)=(1/A(i,i))*(B(i)-sum1-sum2); elseif method==3 X(i)=(OMEGA/A(i,i))*(B(i)-sum1-sum2)-(OMEGA-1)*X(i); end; end;   if method==1 for i=1:n X(i)=Z(i); end end;   precision_flag=0; for i=1:n if (X(i)-previous_X(i) <= TOL) precision_flag=precision_flag+1; end end if (precision_flag==n) break; end     end   if method==1 disp(''); disp('The JACOBI solution is X = '); X disp('Iterations needed : '); disp(repeats); elseif method==2 disp(''); disp('The GAUSS -SEIDEL solution is X = '); X disp('Iterations needed : '); disp(repeats); elseif method==3 disp(''); disp('The SOR solution is X = '); X disp('Iterations needed : '); disp(repeats); end     end   % Under the musical influence of: % "The man who sold the world " (David Bowie)               % % ITERATIVE METHODS % % JACOBI METHOD FOR SOLVING LINEAR NXN SYSTEMS % % % <-Dimitris B. Kalamaras ->   clc; clear all;     disp( 'LINEAR SYSTEMS SOLVE METHODS (iterative)'); disp(' JACOBI method for a NxN linear system ');     n=input( 'ENTER system dimension N:');   disp( 'ENTER matrix A, element by element in a row by row fashion:');   for i=1:n for j=1:n A(i,j)=input ('element : '); end end   disp( 'ENTER matrix B:'); for i=1:n B(i)=input('B element : '); end   disp( 'ENTER matrix X:'); for i=1:n X(i)=input('X element : '); end   max_repeats=input('ENTER maximum repeats : ');     disp( 'MATRIX A:'); A disp('MATRIX B:'); B   %main algorithm for repeats=1:max_repeats for i=1:n sum1=0; sum2=0; if i~=1 for j=1:i-1 sum1=sum1+A(i,j)*X(j); end end   if i~=n for j=i+1:n sum2=sum2+A(i,j)*X(j); end end   Z(i)=(1/A(i,i))*(B(i)-sum1-sum2); end   for i=1:n X(i)=Z(i); end   end       disp(''); disp('The solution is X = '); X       % % ITERATIVE METHODS % % GAUSS SEIDEL METHOD FOR SOLVING LINEAR NXN SYSTEMS % % % <-Dimitris B. Kalamaras -> %     clc; clear all;     disp( 'LINEAR SYSTEMS SOLVE METHODS (iterative)'); disp(' GAUSS SEIDEL method for a NxN linear system ');   n=input( 'ENTER system dimension N:');   disp( 'ENTER matrix A, element by element in a row by row fashion:');   for i=1:n for j=1:n A(i,j)=input ('element : '); end end   disp( 'ENTER matrix B:'); for i=1:n B(i)=input('B element : '); end   disp( 'ENTER matrix X:'); for i=1:n X(i)=input('X element : '); end max_repeats=input('ENTER maximum repeats : ');     disp( 'MATRIX A:'); A disp('MATRIX B:'); B   %main algorithm for repeats=1:max_repeats   for i=1:n sum1=0; sum2=0; if i ~= 1 for j=1:i-1 sum1=sum1+A(i,j)*X(j); end end if i ~= n for j=i+1:n sum2=sum2+A(i,j)*X(j); end end X(i)=(1/A(i,i))*(B(i)-sum1-sum2); end   end       disp(''); disp('The solution is X = '); X           % % ITERATIVE METHODS % % SUCCESIVE OVERRELAXATION METHOD FOR SOLVING LINEAR NXN SYSTEMS % % BEST OMEGA=1.062 :-) % % <-Dimitris B. Kalamaras ->   clc; clear all;     disp( 'LINEAR SYSTEMS SOLVE METHODS (iterative)'); disp(' Succesive Over Relaxation method for a NxN linear system ');   n=input( 'ENTER system dimension N:');   disp( 'ENTER matrix A, element by element in a row by row fashion:');   for i=1:n for j=1:n A(i,j)=input ('element : '); end end   disp( 'ENTER matrix B:'); for i=1:n B(i)=input('B element : '); end   disp( 'ENTER matrix X:'); for i=1:n X(i)=input('X element : '); end max_repeats=input('ENTER maximum repeats : '); OMEGA=input('ENTER relaxation parameter omega : ');     disp( 'MATRIX A:'); A disp('MATRIX B:'); B   %main algorithm for repeats=1:max_repeats   for i=1:n sum1=0; sum2=0; if i ~= 1 for j=1:i-1 sum1=sum1+A(i,j)*X(j); end end if i ~= n for j=i+1:n sum2=sum2+A(i,j)*X(j); end end X(i)=(OMEGA/A(i,i))*(B(i)-sum1-sum2)-(OMEGA-1)*X(i); end   end     disp(''); disp('The solution is X = '); X           % steepest descent with armijo -- Check it . I do not remember if it worked out!!! % % <- Dimitris B. Kalamaras ->   clc; clear all; format long; x(1,1)=input('Enter initial x1(0) [1] : '); x(2,1)=input('Enter initial x2(0) [1]: '); max_repeats=input('Enter maximum repeats [ > 1500 ]:'); tol=input('Enter precision level tol [ < 0.00001 ]: ');     h0=0.99; % this has to be less than 1   for m=1:45 % play with this if you like but 50 is the maximum matlab knows...heheheh hta(m)=h0*2^(1-m); end   h=h0;   for i=2:max_repeats   % anadelta_f=[2-2*x(1,i-1) ; 2-2*x(2,i-1)];   anadelta_f = [2*x(1,i-1) ; 2*x(2,i-1)]; if anadelta_f==0 disp('zero gradient. possible minimum' ); break; end   x(1,i) = x(1,i-1) - h*anadelta_f(1); x(2,i) = x(2,i-1) - h*anadelta_f(2);   f_old = x(1,i-1)^2+x(2,i-1)^2; f_new = x(1,i)^2+x(2,i)^2; %f_old = 2+ 2*x(1,i-1)+ 2*x(2,i-1)- x(1,i-1)^2 - x(2,i-1)^2; %f_new = 2+ 2*x(1,i)+ 2*x(2,i)- x(1,i)^2 - x(2,i)^2;     for m=2:10 if (f_new-f_old) <= -hta(m)/2 * norm(anadelta_f,2)^2 h=hta(m); end end   if abs(x(1,i)-x(1,i-1))             % % LOTKA VOLTERA 2 POPULATION MODEL % % y'=f(x,y,z) = ay-yz % z'=g(x,y,z) = -bz+yz % y0=y(x0) % z0=z(x0) % % Using EULER, EULER-CAUCHY, IMPROVED-EULER, RUNGE-KUTTA 4th CLASS numerical solution methods. % % % if you want to solve another system of ODEs change the f_ode % % <- Dimitris B. Kalamaras ->     clc; clear all;   disp('LOTKA VOLTERA 2 POPULATION MODEL');   disp(''); x(1)=input('Enter x0 [0]: '); y(1)=input('Enter y0 [100]: '); z(1)=input('Enter z0 [15]: '); h=input ('Enter h [0.0005]: '); n=input ('Enter N [3000]: '); a=input ('Enter a [60]: '); b=input ('Enter b [34]: ');   disp('Select a numerical solution ODE method:'); disp('1. Euler'); disp('2. Euler_Cauchy'); disp('3. Improved Euler'); disp('4. Runge 4th Class'); choice=input('> ');     switch (choice) case 1, for i=1:n y(i+1)=y(i)+h*f_ode(x(i),y(i),z(i), 1, a,b); z(i+1)=z(i)+h*f_ode(x(i),y(i),z(i), 2,a,b); x(i+1)=x(i)+h; sprintf('x = %15.10f, y = %15.10f, z = %15.10f',x(i+1),y(i+1),z(i+1)); end hold on; plot(x,y,'r*-'); plot(x,z,'b.:'); hold off; figure; plot(y,z,'g:'); title('LOTKA VOLTERA 2 populations MODEL with EULER method');   case 2,   for i=1:n y(i+1)=y(i)+h*f_ode( x(i)+h/2, y(i)+(h/2)*f_ode(x(i),y(i),z(i),1,a,b), z(i)+(h/2)*f_ode(x(i),y(i),z(i),2,a,b), 1, a,b ); z(i+1)=z(i)+h*f_ode( x(i)+h/2, y(i)+(h/2)*f_ode(x(i),y(i),z(i),1,a,b), z(i)+(h/2)*f_ode(x(i),y(i),z(i),2,a,b), 2, a,b ); x(i+1)=x(i)+h; sprintf('x = %15.10f, y = %15.10f, z = %15.10f',x(i+1),y(i+1),z(i+1));   end hold on; plot(x,y,'r*-'); plot(x,z,'b.:'); hold off; figure; plot(y,z,'g:'); title('LOTKA VOLTERA 2 populations MODEL with EULER-CAUCHY method');   case 3, for i=1:n y(i+1)=y(i)+ (h/2) * ( f_ode( x(i),y(i),z(i),1,a,b ) + f_ode( x(i)+h, y(i)+h*f_ode(x(i),y(i),z(i),1,a,b), z(i)+h*f_ode(x(i),y(i),z(i),2,a,b), 1 , a,b ) ); z(i+1)=z(i)+ (h/2) * ( f_ode( x(i),y(i),z(i),2,a,b ) + f_ode( x(i)+h, y(i)+h*f_ode(x(i),y(i),z(i),1,a,b), z(i)+h*f_ode(x(i),y(i),z(i),2,a,b), 2 , a ,b ) ); x(i+1)=x(i)+h; sprintf('x = %15.10f, y = %15.10f, z = %15.10f',x(i+1),y(i+1),z(i+1)); end hold on; plot(x,y,'r*-'); plot(x,z,'b.:'); hold off; figure; plot(y,z,'g:'); title('LOTKA VOLTERA 2 populations MODEL with IMPROVED EULER');   case 4, for i=1:n k1=h*f_ode(x(i),y(i),z(i),1,a,b); l1=h*f_ode(x(i),y(i),z(i),2,a,b); k2=h*f_ode(x(i)+h/2, y(i)+k1/2, z(i)+l1/2,1,a,b); l2=h*f_ode(x(i)+h/2, y(i)+k1/2, z(i)+l1/2,2,a,b); k3=h*f_ode(x(i)+h/2, y(i)+k2/2, z(i)+l2/2,1,a,b); l3=h*f_ode(x(i)+h/2, y(i)+k2/2, z(i)+l2/2,2,a,b); k4=h*f_ode(x(i)+h, y(i)+k3, z(i)+l3,1,a,b); l4=h*f_ode(x(i)+h, y(i)+k3, z(i)+l3,2,a,b);   y(i+1)=y(i)+ (1/6)*(k1 + 2*k2 + 2*k3 + k4); z(i+1)=z(i)+ (1/6)*(l1 + 2*l2 + 2*l3 + l4);   x(i+1)=x(i)+h; sprintf('x = %15.10f, y = %15.10f, z = %15.10f',x(i+1),y(i+1),z(i+1)); %exclamation mark here to abandon output. end hold on; plot(x,y,'r*-'); plot(x,z,'b.:'); hold off; figure; plot(y,z,'g:'); title('LOTKA VOLTERA 2 populations MODEL with RUNGE KUTTA 4TH CLASS');   otherwise for i=1:n k1=h*f_ode(x(i),y(i),z(i),1,a,b); l1=h*f_ode(x(i),y(i),z(i),2,a,b); k2=h*f_ode(x(i)+h/2, y(i)+k1/2, z(i)+l1/2,1,a,b); l2=h*f_ode(x(i)+h/2, y(i)+k1/2, z(i)+l1/2,2,a,b); k3=h*f_ode(x(i)+h/2, y(i)+k2/2, z(i)+l2/2,1,a,b); l3=h*f_ode(x(i)+h/2, y(i)+k2/2, z(i)+l2/2,2,a,b); k4=h*f_ode(x(i)+h, y(i)+k3, z(i)+l3,1,a,b); l4=h*f_ode(x(i)+h, y(i)+k3, z(i)+l3,2,a,b);   y(i+1)=y(i)+ (1/6)*(k1 + 2*k2 + 2*k3 + k4); z(i+1)=z(i)+ (1/6)*(l1 + 2*l2 + 2*l3 + l4);   x(i+1)=x(i)+h; sprintf('x = %15.10f, y = %15.10f, z = %15.10f',x(i+1),y(i+1),z(i+1)); %exclamation mark here to abandon output. end hold on; plot(x,y,'r*-'); plot(x,z,'b.:'); hold off; figure; plot(y,z,'g:'); title('LOTKA VOLTERA 2 populations MODEL with RK 4TH CLASS ');   end       function y=f_ode(x1,x2,x3,choice,a,b)   % function 'f_ode' is the the functions % used in the LOTKA VOLTERA MODEL with RK, EULER, EULER-CAUCHY (above) % % <- Kalamaras Dimitris ->     if nargin==0 error('Usage: f_ode(x1,x2,x3,choice,a,b ) , where xi real, choice {1 or 2} and a,b parameters'); end   if choice==1 y=a*x2 - x2*x3 ; %% lotka VOLTERA MODEL f(x,y,z)=ay-yz %% elseif choice==2 y=-b*x3+ x2*x3; %% lotka voltera MODEL g(x,y,z)=-bz+yz %% else error('Usage: f_ode(x1,x2,x3,choice,a,b ) , where xi real, choice {1 or 2} and a,b parameters'); end         % Approximation of the integral of 4/(1+x^2) from zero to one via Runge - Kutta 4th class % % set y'=4/(1+x^2) and history y(0)=0, thus the integral is the solution of this ODE in 0<=x<=1 if y(1)=pi % where pi is the real value of the given integral. The method used in the integration of the ODE is Runge Kutta % % <- Kalamaras Dimitris ->     clc; clear all; format long;   h=input ('Enter step h for R-K : '); %h=0.01; tol=input('Enter desired precision TOL : ' );   max=1/h; x=0; % y=0; % history zero; z=0; for i=1:max   k1=h*f_rk(x,y); k2=h*f_rk(x+h/2,y+k1/2); k3=h*f_rk(x+h/2,y+k2/2); k4=h*f_rk(x+h,y+k3); y=y+ (1/6) * (k1+2*k2+2*k3+k4); if abs(y-z)<= tol break; end z=y; if x == 1 disp('0k');break; end x=i*h;     end disp(' '); disp('RUNGE KUTTA RESULT:'); sprintf(' y(1) = %15.14f in %5d RUNGE KUTTA iterations', y, i-1)   disp('while the real integral(4/(1+x^2)) is ' ); syms r; real_value=int(4/(1+r^2),0,1)     difference=abs(real_value-y)     x=0; y=0; z=0; for i=1:max   k1=h*f_rk(x,y); k2=h*f_rk(x+h,y+k1); y=y+ (1/2)*(k1+k2); if abs(y-z)<= tol break; end z=y; if x == 1 disp('0k');break; end x=i*h;     end   disp(' ' ); disp('EULER - HEUN RESULT:'); sprintf(' y(1) = %15.14f in %5d EULER-HEUN iterations', y, i-1)   difference=abs(y-real_value)       function z=f_rk(x,y)   % function 'f_rk' is the the function % used in the RungeKutta method above % z=4/(1+x^2);             % Boundary Value D.E % Numerical Solution usinh the Shooting non-linear method with RUNGE KUTTA 4th class for % the solution of the initial value problems and Bisection and Newton Raphson method for t % % y''=f(x,y,y') = 3y'-2y+x % if you want another system, change f_ode1.m % y(0)=1 % y(1)=2 shooting: y'(0)=t % % % <- Dimitris B. Kalamaras ->   clc; clear all;   disp( 'SHOOTING METHOD with RUNGE-KUTTA 4th Class and NEWTON-RAPHSON or BISECTION '); disp( ' y''''=f(x,y,y'') = 3y''-2y+x'); disp( ' y(0)=1, y(1)=1 --> shooting: y''(0)=t'); disp( ' ');   %DATA INPUT & INITIALIZATION h=input ('RUNGE-KUTTA...: Enter step h for RK [0.001]: '); e=input ('NEWTON-RAPHSON: Enter step e for NR [0.1]: '); t(1)=input('NR: TARGET....: Enter initial estimation of t : '); % newton raphson only prec =input('NR & BISECTION: t approximation precision [<0.1]: '); a=input ('BISECTION.... : Enter left interval margin a [0.1]: '); b=input ('BISECTION ....: Enter right interval margin b [3]: '); max=input ('Enter max. amount of iterations for shooting [8] : ') ; n=1/h; Target=2;     % % Runge Kutta with Finite-Differences Newton-Raphson % for k=1:max % solve the equation P(t)=y(b;t)-B=0 with a method, f.ex. Newton Raphson to find a better t. % Thus we need the root of the P % I need the derivative Of P(t)=y(b,t)-B for Newton Raphson, thus i use dP/dt ~ ( P(t+e)-P(t) )/e as an approxim. % for that reason we calculate first y with t=t+e: x(1)=0; %bored changing the damn zero-index matlab option...who cares whilst it works, ? y(1)=1; z(1)=t(k)+e;   for i=1:n k1=h*f_ode1(x(i),y(i),z(i),1); l1=h*f_ode1(x(i),y(i),z(i),2); k2=h*f_ode1(x(i)+h/2, y(i)+k1/2, z(i)+l1/2,1); l2=h*f_ode1(x(i)+h/2, y(i)+k1/2, z(i)+l1/2,2); k3=h*f_ode1(x(i)+h/2, y(i)+k2/2, z(i)+l2/2,1); l3=h*f_ode1(x(i)+h/2, y(i)+k2/2, z(i)+l2/2,2); k4=h*f_ode1(x(i)+h, y(i)+k3, z(i)+l3,1); l4=h*f_ode1(x(i)+h, y(i)+k3, z(i)+l3,2);   y(i+1)=y(i)+ (1/6)*(k1 + 2*k2 + 2*k3 + k4); z(i+1)=z(i)+ (1/6)*(l1 + 2*l2 + 2*l3 + l4);   x(i+1)=x(i)+h;   end   Pte=(y(n+1)-Target);   % Now calculate y with target t=t   x(1)=0; %bored changing the damn zero-index matlab option...who cares whilst it works, ? y(1)=1; z(1)=t(k);   for i=1:n k1=h*f_ode1(x(i),y(i),z(i),1); l1=h*f_ode1(x(i),y(i),z(i),2); k2=h*f_ode1(x(i)+h/2, y(i)+k1/2, z(i)+l1/2,1); l2=h*f_ode1(x(i)+h/2, y(i)+k1/2, z(i)+l1/2,2); k3=h*f_ode1(x(i)+h/2, y(i)+k2/2, z(i)+l2/2,1); l3=h*f_ode1(x(i)+h/2, y(i)+k2/2, z(i)+l2/2,2); k4=h*f_ode1(x(i)+h, y(i)+k3, z(i)+l3,1); l4=h*f_ode1(x(i)+h, y(i)+k3, z(i)+l3,2);   y(i+1)=y(i)+ (1/6)*(k1 + 2*k2 + 2*k3 + k4); z(i+1)=z(i)+ (1/6)*(l1 + 2*l2 + 2*l3 + l4);   x(i+1)=x(i)+h; end disp('SHOOTING with NEWTON-RAPHSON:'); sprintf('In %2d repeat, for n= %6d the desired y(n,t) has value %20.17f with t = %8.4f', k, n, y(n+1),t(k))   Pt=y(n+1)-Target; % function P(t) of the difference y(b,t)-target will be used in Newton Raphson to find a better t if abs(Pt) < prec disp('Shooting Succeeded!'); t(k+1)=t(k); break; end;   % approximation of the derivative dPt=(Pte-Pt)/e; t(k+1)=t(k)-Pt/dPt;   figure; plot(x,y,'g:'); title('SHOOTING & NR: Graph of y in between approximations of t '); end figure; plot(x,y,'r:'); title('SHOOTING & NR: Final graph of y in the interval [0,1]'); disp('NR Iterations Needed :' ); disp (k);       disp('Wait...Writing plot for Newton-Raphson...'); disp(' ');               % % Runge Kutta with BISECTION %   % FIRST CALCULATE y with target t=a , t=b, i.e. the left % right margins of the interval... % this is done because I need the value (sign) of the function Pt= y(b,t)-Target at t = a % and also the value at t=b for k=1:2   x(1)=0; y(1)=1; if k==1 z(1)=a; elseif k==2 z(1)=b; end   for i=1:n k1=h*f_ode1(x(i),y(i),z(i),1); l1=h*f_ode1(x(i),y(i),z(i),2); k2=h*f_ode1(x(i)+h/2, y(i)+k1/2, z(i)+l1/2,1); l2=h*f_ode1(x(i)+h/2, y(i)+k1/2, z(i)+l1/2,2); k3=h*f_ode1(x(i)+h/2, y(i)+k2/2, z(i)+l2/2,1); l3=h*f_ode1(x(i)+h/2, y(i)+k2/2, z(i)+l2/2,2); k4=h*f_ode1(x(i)+h, y(i)+k3, z(i)+l3,1); l4=h*f_ode1(x(i)+h, y(i)+k3, z(i)+l3,2); y(i+1)=y(i)+ (1/6)*(k1 + 2*k2 + 2*k3 + k4); z(i+1)=z(i)+ (1/6)*(l1 + 2*l2 + 2*l3 + l4);   x(i+1)=x(i)+h; end   if k==1 Pta=(y(n+1)-Target); % Got it. elseif k==2 Ptb=(y(n+1)-Target); % Got it. end;   end   t(1)=a+((b-a)/2^1)*sign(Pta)*sign(Pta); % calculate initial t with BISECTION   % Now approximate t(k) with Bisection in max times, by calculating y, and testing the Pt=y(b,t)-Target function.   for k=1:max   % Now calculate y with target t=t(k) given from Bisection   x(1)=0; %bored changing the damn zero-index matlab option...who cares whilst it works, ? y(1)=1; z(1)=t(k);   for i=1:n k1=h*f_ode1(x(i),y(i),z(i),1); l1=h*f_ode1(x(i),y(i),z(i),2); k2=h*f_ode1(x(i)+h/2, y(i)+k1/2, z(i)+l1/2,1); l2=h*f_ode1(x(i)+h/2, y(i)+k1/2, z(i)+l1/2,2); k3=h*f_ode1(x(i)+h/2, y(i)+k2/2, z(i)+l2/2,1); l3=h*f_ode1(x(i)+h/2, y(i)+k2/2, z(i)+l2/2,2); k4=h*f_ode1(x(i)+h, y(i)+k3, z(i)+l3,1); l4=h*f_ode1(x(i)+h, y(i)+k3, z(i)+l3,2);   y(i+1)=y(i)+ (1/6)*(k1 + 2*k2 + 2*k3 + k4); z(i+1)=z(i)+ (1/6)*(l1 + 2*l2 + 2*l3 + l4);   x(i+1)=x(i)+h; end disp('Shooting with BISECTION :'); sprintf('In %2d repeat, for n= %6d the desired y(n,t) has value %20.17f with t = %8.4f', k, n, y(n+1),t(k))   Ptk=y(n+1)-Target; % function P(t(k)) of their difference will be used in BISECTION to find a better t if abs(Ptk) < prec disp('Shooting Succeeded!'); t(k+1)=t(k); break; end;     % BISECTION ALGORITHM   if Pta*Ptk >0 & Ptb*Ptk > 0 ... disp('No roots for P(t)=y(b,t)-Target in the given interval '); break; end;   t(k+1)=t(k)+((b-a)/2^(k+1))*sign(Pta)*sign(Ptk);       figure; plot(x,y,'g:'); title('SHOOTING & BISECTION: Graph of y in between approximations of t '); end figure; plot(x,y,'r:'); title('SHOOTING & BISECTION: Final graph of y in the interval [0,1]'); disp('BISECTION Iterations Needed :' ); disp (k);         %for educational purposes plot y from 0 to 5   x(1)=0; y(1)=1; z(1)=t(k+1); h=0.001; n=5/h ; for i=1:n k1=h*f_ode1(x(i),y(i),z(i),1); l1=h*f_ode1(x(i),y(i),z(i),2); k2=h*f_ode1(x(i)+h/2, y(i)+k1/2, z(i)+l1/2,1); l2=h*f_ode1(x(i)+h/2, y(i)+k1/2, z(i)+l1/2,2); k3=h*f_ode1(x(i)+h/2, y(i)+k2/2, z(i)+l2/2,1); l3=h*f_ode1(x(i)+h/2, y(i)+k2/2, z(i)+l2/2,2); k4=h*f_ode1(x(i)+h, y(i)+k3, z(i)+l3,1); l4=h*f_ode1(x(i)+h, y(i)+k3, z(i)+l3,2);   y(i+1)=y(i)+ (1/6)*(k1 + 2*k2 + 2*k3 + k4); z(i+1)=z(i)+ (1/6)*(l1 + 2*l2 + 2*l3 + l4);   x(i+1)=x(i)+h;   end figure; plot(x,y); title('BISECTION Final Plot of integral y of the given ODE for the interval [0,5] ');       function y=f_ode1(x1,x2,x3,choice)   % function 'f_ode1' is the system of the functions F,G % used in the with RK above % % The initial valued ODE { y''=3y'-2y+x y(0)= 1, y'(0)=t } % with y1=y, y2=y', y3=y'' becomes the system of 1st class ODEs: % -- { y1'=y2, with y1(0)=1 % -- { y2'=y''=3y'-2y+x=3y2-2y+x with y2(0)=t % % <- Kalamaras Dimitris ->         if nargin==0 error('Usage: f_ode1(x1,x2,x3,choice ) , where xi real, choice {1 or 2} '); end   if choice==1 y=x3 ; %% %% elseif choice==2 y=3*x3-2*x2+x1; %% %% else error('Usage: f_ode1(x1,x2,x3,choice ) , where xi real, choice {1 or 2} '); end         * RANDOM NUMBER GENERATOR ADAPTED FROM NUMERICAL RECIPES   function y=ran0(seed)   % Returns a uniform random deviate between 0.0 and 1.0. % The period of ran0 is 2^31~2.1x10^9 % % Set or reset idum to any integer (except the unlikely value MASK) to initialize % the sequence; idum must not be altered between calls for successive deviates % in a sequence. Value 0 must not be allowed as seed and should not be expected as rand.num! % This is a peculiarity of all r.n.g of the form I(n+1)=a*I(n) (mod m) % % Taken from Numerical Recipes in C Internet Edition. % Adapted in Matlab for creating noise over a function. % View ran0.m for more details...if you have lost NR in C! % % <- Kalamaras Dimitris ->   format long; global idum;   % Park & Miller "Minimal Standard" proposal if nargin==1 idum=seed; end idum; ia=16807; im=2147483647; %2^31-1 am=1.0/im; %MASK=123459876 % if im=ia*iq+ir (Shcrage Method) % then ia*idum=a(z mod q)-r[z/q] (+im if <0) % for multiplying two 32-bit integers modulo a 32-bit constant, % without using any intermediates larger than 32 bits (including a sign bit)   iq=127773; ir=2836; % I should XOR with MASK to allow zero seed (read NR) but... %idum=xor(idum,MASK) % unfortunately XORing in matlab requires vectors x,y! k=round(idum/iq); % Done with round (thus, less efficiently) to keep code NR-like idum=ia*(idum-iq*k)-ir*k; % It also seems that I am avoidind overflows...DO NOT USE MOD! if idum < 0 % If it is <0 im should be added. idum=idum+im; end %convert to fp y=am*idum; %idum=xor(idum,MASK)     *   % NUMERICAL ANALYSIS % This script integrates error function of a standard normal variable % manifestated in f_mc_erf.m using the Mean Sample Monte Carlo method % % <- Kalamaras Dimitris ->     clc; clear all; format long; disp('NUMERICAL INTEGRATION'); disp('Monte Carlo Method of error function in (0,1)'); tic; I_Estimator=zeros(5,1); %n=input('Enter Amount of Random Draws for Monte Carlo :'); n=40000; trial=1; while n < 40001 sum=0; rand('state',0); for i=1:n x=rand; sum=sum+F_mc_erf(x); end   I_Estimator(trial)=(1/n)*sum; sum=0; rand('state',0); for i=1:n x=rand; sum=sum+(F_mc_erf(x)-I_Estimator(trial))^2; end   disp(sum);   Var_Estimator=sqrt( (1/n) * (1/(n-1))*sum ); Var_1=sqrt ( (1/n) * (1/(n-1)) * sum ); %??????????????????????? zero ????????   disp(n); sprintf('The approximate value of the integral is I = %10.5f',I_Estimator(trial)) sprintf('With Error of the I to be = %10.5f',Var_Estimator) sprintf('With Variance of the I to be = %10.5f',Var_1) n=n+1000; trial=trial+1; end     x=zeros(5,1);   for i=1:5 if i==1 x(i)=1000 else x(i)=x(i-1)+1000; end end   plot (x, I_Estimator, '.'); toc;               % MONTE CARLO INTEGRATION % This script integrates a given univariate function manifestated in f_mc1.m % using the Monte Carlo method % % <- Kalamaras Dimitris ->     clc; clear all; disp('NUMERICAL INTEGRATION'); disp('Monte Carlo Method'); n=input('Enter Amount of Random Draws for Monte Carlo :'); x=zeros(n); sum=0;   for i=1:n x(i)=rand; sum=sum+F_mc1(x(i)); end   I_Estimator=(1/n)*sum; sum=0; for i=1:n sum=sum+(F_mc1(x(i))-I_Estimator)^2; end   Var_Estimator=(1/n)*(1/(n-1))*sum;   sprintf('The approximate value of the integral is I = %10.5f',I_Estimator) sprintf('With Variance of the I to be = %10.5f',Var_Estimator)       function y=f_mc1(x)   % function 'f_mc1' is the the host m file of the function which % is being integrated by the "mc1" script % % <- Kalamaras Dimitris ->     if nargin==0 x=0:.1:10; end %if nargin==0 error('Usage: f_mc1 ( x ) , where x whatever.'); end y=cos(x); %y=x+x^4;                 % This script integrates a given two-variate function manifestated in f_mc2.m % using the Monte Carlo method   % <- Kalamaras Dimitris ->     clc; clear all; disp('NUMERICAL INTEGRATION'); disp('Monte Carlo Method in 2 dimensions '); a1=input('Input left margin a1 for the first dimension :'); b1=input('Input right margin b1 for the first dimension :'); a2=input('Input left margin a2 for the second dimension :'); b2=input('Input right margin b2 for the second dimension :'); k=1/((b1-a1)*(b2-a2)) n=input('Enter Amount of Random Draws for Monte Carlo :'); %x=zeros(n); %y=zeros(n); sum=0;   for i=1:n % x(i)=(b1-a1)*rand-a1; %transformers....enotheite! % y(i)=(b2-a2)*rand-a2; x=(b1-a1)*rand-a1; y=(b2-a2)*rand-a2; % sum=sum+F_mc2(x(i),y(i)); sum=sum+F_mc2(x,y);   end   I_Estimator=(1/(k*n))*sum;   %sum=0; %for i=1:n % sum=sum+(F_mc1(x(i))-I_Estimator)^2; %end   %Var_Estimator=(1/n)*(1/(n-1))*sum;   sprintf('The approximate value of the integral is I = %10.5f',I_Estimator) %sprintf('With Variance of the I to be = %10.5f',Var_Estimator)       function y=f_mc2(x,y)   % function 'f_mc2' is the the host m file of the 2dimensional function which % is being integrated by the "mc2" script % % <- Kalamaras Dimitris ->       if (nargin==0) x=0:.1:10; y=0:.1:10; end %if nargin==0 error('Usage: f_mc2 ( x,y ) , where x whatever.'); end     %y=x^2+y^2*cos(x)-exp(x); y=x^2+y^2;```