Analysis
September 27, 2008
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;
