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; |
Leave a Reply