MatLab homework - root solving
MATLAB/Assignment 2/Bracketing (Multiple Roots) (4)/Bisection method bracketing/Thumbs.db
MATLAB/Assignment 2/Bracketing (Multiple Roots) (4)/Newton method bracketing/Thumbs.db
MATLAB/Assignment 2/Bracketing (Multiple Roots) (4)/Thumbs.db
MATLAB/Assignment 2/Thumbs.db
MATLAB/Assignment 2/Newton method - intersection of two functions (3)/Thumbs.db
MATLAB/Assignment 2/Bisecting method (1)/Thumbs.db
MATLAB/Assignment 2/Newton Method (2)/Thumbs.db
MATLAB/Assignment 1/Images/Thumbs.db
MATLAB/Assignment 1/Nano Sim Assignment 1.docx
Simulation on the Nanoscale ENG-M03: Assignment 1
Using a script file (I)
clear
format short g
N = 100;
x = linspace(-pi,pi,N)';
y = sin(x);
figure(1), plot(x,y)
Defining a function (II)
function y = sin_scale (x,sf)
y = sin(x./sf);
end
sf = 0.1;
y = sin_scale(x,sf);
figure(2), plot(x,y)
sf=[0.1 0.5 2];
N_sf = numel(sf);
sin_array = zeros(N,N_sf);
for ndx = 1: N_sf
sin_array(:,ndx) = sin_scale(x,sf(ndx));
end
figure(3),plot(x,sin_array)
Numerical Differentiation
Function
function dydx = Deriv(x,y)
N = numel(x);
dydx = zeros(N-1,1);
h= x(2)-x(1);
for ndx = 1 : N-1
dydx(ndx,:) =(y(ndx+1)-y(ndx))/h;
end
Script
clear
format short g
N=100;
x = linspace (-pi,pi,N)';
y = cos(x);
dydx = Deriv(x,y);
xd = x(1:N-1);
figure(4), plot(x,y,xd,dydx,x,-sin(x))
Numerical Integration
1. Trapezoidal Method
a) Integral
b) g=integral(@f,a,b,'ArrayValued',true)
Error = ((I-g)/g).*100 = 23.37%
2. Composite Trapezium Rule
a) function [ aproxInt ] = aproxInt(n,a,b)
function y = sumFx
y = 0;
for k = 1:n-1
y = y + feval(@f,(a+(k.*((b-a)/n))));
end
end
Totx=sumFx;
aproxInt = 0.5.*((b-a)/n).*((feval(@f,a))+(feval(@f,b))+(2.*(Totx)));
end
b) I= aproxInt(n,a,b)= 1.00002056192951
Error = ((I-g)/g).*100 = 0.00205619295078341
c) I1 = aproxInt2(5,1,2) = 0.695634920634921
I2 = aproxInt2(50,1,2) = 0.693172179310195
I3 = aproxInt2(100,1,2) = 0.693153430481824
d) g1 = log(2);
e1 = ((I1-g1)/g1).*100 = 0.358905026918763
e2 = ((I2-g1)/g1).*100 = 0.00360655730140557
e3 = ((I3-g1)/g1).*100 = 0.000901673130050152
e) number of sub-intervals needed to achieve an accuracy to 6 d.p., n=277 (adjusted manually)
3. Simpson’s Rule
MATLAB/Assignment 2/Nanosim assignment 2.docx
Simulation on the Nanoscale ENG-M03: Assignment 2
1. Bisection Method
2. Newton Method
3. Newton Method – Intersection of Two Functions
This script only takes 3 iterations to calculate the root. Newton’s method by intersection of two functions is much faster than the bisection method, but encounters problems unless the root is well defined (e.g. the wrong root may be found).
4. Bracketing (Bisection Method)
MATLAB/Assignment 1/Composite Trapezium Rule(5)/aproxInt.m
%5a) function [ aproxInt ] = aproxInt(n,a,b) function y = sumFx y = 0; for k = 1:n-1 y = y + feval(@f,(a+(k.*((b-a)/n)))); end end Totx=sumFx; aproxInt = 0.5.*((b-a)/n).*((feval(@f,a))+(feval(@f,b))+(2.*(Totx))); end
MATLAB/Assignment 1/Composite Trapezium Rule(5)/aproxInt2.m
function [ aproxInt2 ] = aproxInt2(n,a,b) function y = sumGx y = 0; for k = 1:n-1 y = y + feval(@g,(a+(k.*((b-a)/n)))); end end Totx=sumGx; aproxInt2 = 0.5.*((b-a)/n).*((feval(@g,a))+(feval(@g,b))+(2.*(Totx))); end
MATLAB/Assignment 2/Bisecting method (1)/bisect.m
function root = bisect(f,a,b,imax,tol) % Input: f function to be evaluated % a,b interval limits % imax iteration maximum % tol convergence criterion % Output: root position of the root % initialise parameters, i.e. x1=a; f1=feval(f,x1); x3=b; f3=feval(f,x3); % CHECK: verify that there is indeed a root in the interval [a b]. if f1*f3>0 error('No root in the specified interval') end for ndx=1:imax x2=(x3+x1)/2; f2=feval(f,x2); %etc… % Find which half interval contains root – i.e. if-else construction % test for convergence i.e. d < tol % CHECK: Excessive iterations if f1*f2<0 % Root is in left half, so d=(x2+x1)/2; f3=f2; x3=x2; else % Root is in right half, so d=(x2+x3)/2; f1=f2; x1=x2; end d=abs(x1-x2); if d<tol root =x2; end end end
MATLAB/Assignment 2/Bracketing (Multiple Roots) (4)/Bisection method bracketing/bisect.m
function root = bisect(f,a,b,imax,tol) % Input: f function to be evaluated % a,b interval limits % imax iteration maximum % tol convergence criterion % Output: root position of the root % initialise parameters, i.e. x1=a; f1=feval(f,x1); x3=b; f3=feval(f,x3); % CHECK: verify that there is indeed a root in the interval [a b]. if f1*f3>0 error('No root in the specified interval') end for ndx=1:imax x2=(x3+x1)/2; f2=feval(f,x2); %etc… % Find which half interval contains root – i.e. if-else construction % test for convergence i.e. d < tol % CHECK: Excessive iterations if f1*f2<0 % Root is in left half, so d=(x2+x1)/2; f3=f2; x3=x2; else % Root is in right half, so d=(x2+x3)/2; f1=f2; x1=x2; end d=abs(x1-x2); if d<tol root =x2; end end end
MATLAB/Assignment 2/Bisecting method (1)/Bisect_script.m
%1b) clear format short g a=3; b=10; imax=50; tol=1e-6; f ='function1'; xr = bisect(f,a,b,imax,tol); x = linspace(a,b,100)'; y = feval(f,x); yr = feval(f,xr); %figure(1);plot(x,y,xr,yr,'o'); %title('X root = 3.1416'); %1c) clear format short g a=4; b=4.7; imax=50; tol=1e-6; f ='function2'; xr=bisect(f,a,b,imax,tol); x = linspace(a,b,100)'; y = feval(f,x); yr = feval(f,xr); figure(2);plot(x,y,xr,yr,'o'); title('X root = 4.2711');
MATLAB/Assignment 2/Bracketing (Multiple Roots) (4)/Bisection method bracketing/Bisect_script.m
clear format short g imax=50; tol=1e-6; f ='function1'; a=3; b=5; xr1 = bisect(f,a,b,imax,tol); a=5.7; b=6.5; xr2 = bisect(f,a,b,imax,tol); a=7.5; b=8.5; xr3 = bisect(f,a,b,imax,tol); a=9.5; b=10.4; xr4 = bisect(f,a,b,imax,tol); x = linspace(4,10,100)'; f ='function2'; yr1 = feval(f,xr1); yr2 = feval(f,xr2); yr3 = feval(f,xr3); yr4 = feval(f,xr4); y = feval(f,x); f ='function3'; y1 = feval(f,x); figure(1);plot(x,y,'-',x,y1,'-',xr1,yr1,'v',xr2,yr2,'v',xr3,yr3,'v',xr4,yr4,'v'); %title('X root = 3.1416'); %1c) %clear %format short g %a=4; %b=4.7; %imax=50; %tol=1e-6; %f ='function2'; %xr=bisect(f,a,b,imax,tol); %x = linspace(a,b,100)'; %y = feval(f,x); %yr = feval(f,xr); %figure(2);plot(x,y,xr,yr,'o'); %title('X root = 4.2711');
MATLAB/Assignment 1/Lab_1 Q 1-4/Deriv.m
function dydx = Deriv(x,y) N = numel(x); dydx = zeros(N-1,1); h= x(2)-x(1); for ndx = 1 : N-1 dydx(ndx,:) =(y(ndx+1)-y(ndx))/h; end
MATLAB/Assignment 1/Simpsons rule(6)/f.m
function [ f ] = f( x ) f=1/x; end
MATLAB/Assignment 1/Composite Trapezium Rule(5)/f.m
function [ f ] = f( x ) f=x.*sin(x); end
MATLAB/Assignment 1/lab1_trap2 (5)/f.m
function [ f ] = f( x ) f=x.*sin(x); end
MATLAB/Assignment 1/Trapezoidal Method(4)/f.m
function f = f(x) %UNTITLED13 Summary of this function goes here % Detailed explanation goes here f=x'*sin(x); end
MATLAB/Assignment 2/Bracketing (Multiple Roots) (4)/Newton method bracketing/f1_deriv.m
function df1 = f1_deriv (x) df1 = -100.*exp(-x) - (5*pi/2)*cos((pi*x)/2); end
MATLAB/Assignment 2/Newton method - intersection of two functions (3)/f1_deriv.m
function df1 = f1_deriv (x) df1 = -100.*exp(-x) - (5*pi/2)*cos((pi*x)/2); end
MATLAB/Assignment 2/Bracketing (Multiple Roots) (4)/Newton method bracketing/f2_deriv.m
function df1 = f2_deriv (x) df1 = -100.*exp(-x); end
MATLAB/Assignment 2/Newton method - intersection of two functions (3)/f2_deriv.m
function df1 = f2_deriv (x) df1 = -100.*exp(-x); end
MATLAB/Assignment 2/Bracketing (Multiple Roots) (4)/Newton method bracketing/f3_deriv.m
function df3 = f3_deriv (x) df3 = (5*pi/2)*cos((pi*x)/2); end
MATLAB/Assignment 2/Newton method - intersection of two functions (3)/f3_deriv.m
function df3 = f3_deriv (x) df3 = (5*pi/2)*cos((pi*x)/2); end
MATLAB/Assignment 2/Newton Method (2)/func1_deriv.m
function df= func1_deriv(x) df = exp(x) - (1 + x); end
MATLAB/Assignment 2/Newton Method (2)/func2_deriv.m
function df2 = func2_deriv(x) df2 = 10/x - 1; end
MATLAB/Assignment 2/Bisecting method (1)/function1.m
function y = function1(x) y = sin(x).*sinh(x); end
MATLAB/Assignment 2/Newton Method (2)/function1.m
function f = function1(x) f = exp(x) - (1 + x + (x.^2)./2); end
MATLAB/Assignment 2/Bracketing (Multiple Roots) (4)/Bisection method bracketing/function1.m
function f = function1(x) f = 100.*exp(-x) - 5.*sin((pi*x)/2); end
MATLAB/Assignment 2/Bracketing (Multiple Roots) (4)/Newton method bracketing/function1.m
function f = function1(x) f = 100.*exp(-x) - 5.*sin((pi*x)/2); end
MATLAB/Assignment 2/Newton method - intersection of two functions (3)/function1.m
function f = function1(x) f = 100.*exp(-x) - 5.*sin((pi*x)/2); end
MATLAB/Assignment 2/Bracketing (Multiple Roots) (4)/Bisection method bracketing/function2.m
function f2 = function2 (x) f2 = 100.*exp(-x); end
MATLAB/Assignment 2/Bracketing (Multiple Roots) (4)/Newton method bracketing/function2.m
function f2 = function2 (x) f2 = 100.*exp(-x); end
MATLAB/Assignment 2/Newton method - intersection of two functions (3)/function2.m
function f2 = function2 (x) f2 = 100.*exp(-x); end
MATLAB/Assignment 2/Newton Method (2)/function2.m
function f2 = function2(x) f2 = 10.*log(x) - x; end
MATLAB/Assignment 2/Bisecting method (1)/function2.m
function f = function2(x) R = 10; f = sqrt(R.^2 - x.^2) - x.*tan(x); end
MATLAB/Assignment 2/Bracketing (Multiple Roots) (4)/Bisection method bracketing/function3.m
function f3 = function3 (x) f3 = 5.*sin((pi*x)/2); end
MATLAB/Assignment 2/Bracketing (Multiple Roots) (4)/Newton method bracketing/function3.m
function f3 = function3 (x) f3 = 5.*sin((pi*x)/2); end
MATLAB/Assignment 2/Newton method - intersection of two functions (3)/function3.m
function f3 = function3 (x) f3 = 5.*sin((pi*x)/2); end
MATLAB/Assignment 1/Composite Trapezium Rule(5)/g.m
function [ g ] = g( x ) g=1/x; end
MATLAB/Assignment 2/Newton method - intersection of two functions (3)/Intersection_Script.m
clear format short g x=4; imax=20; tol=1e-10; f ='function1'; df ='f1_deriv'; xr = Newton(x,f,df,imax,tol); x = linspace(3,10,100)'; f ='function2'; yr = feval(f,xr); y1 = feval(f,x); f ='function3'; y2 = feval(f,x); yr = feval(f,xr); figure(1);plot(x,y1,'-',x,y2,'-',xr,yr,'o'); title('x root = 4.1949')
MATLAB/Assignment 2/Bracketing (Multiple Roots) (4)/Newton method bracketing/Intersection_Script.m
clear format short g imax=20; tol=1e-6; f ='function1'; df ='f1_deriv'; %roots in the interval x [4 10], approx at [4,1.5][6,0][8,0][10,0] x=4; xr1 = Newton(x,f,df,imax,tol); x=6; xr2 = Newton(x,f,df,imax,tol); x=8; xr3 = Newton(x,f,df,imax,tol); x=10; xr4 = Newton(x,f,df,imax,tol); x = linspace(4,10,100)'; f ='function2'; yr1 = feval(f,xr1); yr2 = feval(f,xr2); yr3 = feval(f,xr3); yr4 = feval(f,xr4); y = feval(f,x); f ='function3'; y1 = feval(f,x); figure(1);plot(x,y,'-',x,y1,'-',xr1,yr1,'o',xr2,yr2,'o',xr3,yr3,'o',xr4,yr4,'o'); %title('x root = 4.1949')
MATLAB/Assignment 1/Lab_1 Q 1-4/lab1_sin.m
clear format short g N = 100; x = linspace(-pi,pi,N)'; y = sin(x); %figure(1), plot(x,y); sf = 0.1; y = sin_scale(x,sf); %figure(2), plot(x,y); sf=[0.1 0.5 2]; N_sf = numel(sf); sin_array = zeros(N,N_sf); for ndx = 1: N_sf sin_array(:,ndx) = sin_scale(x,sf(ndx)); end %figure(3),plot(x,sin_array) %title(['A Sine fn; limits = [-pi pi]; Scale factors = ',num2str(sf)]) clear format short g N=100; x = linspace (-pi,pi,N)'; y = cos(x); dydx = Deriv(x,y); xd = x(1:N-1); %figure(4), plot(x,y,xd,dydx,x,-sin(x)) clear format short g N = 100; x = linspace(-pi,pi,N)'; a = 0; b = pi/2; h = (b-a); f = x.*sin(x); int = trap1(f,a,b);
MATLAB/Assignment 1/Trapezoidal Method(4)/lab1_trap1.m
clear format short g N=100; x=linspace(-pi,pi,N); a=0; b=pi/2; I=trap1(@f,a,b) g=integral(@f,a,b,'ArrayValued',true) Error = ((I-g)/g).*100
MATLAB/Assignment 1/lab1_trap2 (5)/lab1_trap2.m
clear format short g N=100; x=linspace(-pi,pi,N); a=0; b=pi/2; I=trap2(@f,a,b) g=integral(@f,a,b,'ArrayValued',true); Error = ((I-g)/g).*100
MATLAB/Assignment 1/Composite Trapezium Rule(5)/lab1_trap2.m
clear format long g b=pi/2; a=0; n=100; % 5b) s=aproxInt(n,a,b); g=integral(@f,a,b,'ArrayValued',true); Error = ((s-g)/g).*100; % 5c) s1 = aproxInt2(5,1,2); s2 = aproxInt2(50,1,2); s3 = aproxInt2(100,1,2); % 5d) g1 = log(2); e1 = ((s1-g1)/g1).*100; e2 = ((s2-g1)/g1).*100; e3 = ((s3-g1)/g1).*100; % 5) I just adjusted the n value manually until i hit the right number s4 = aproxInt2(277,1,2);
MATLAB/Assignment 1/Simpsons rule(6)/lab1_trap3.m
clear format long g % 6b) s1 = simpsonInt(5,1,2) s2 = simpsonInt(50,1,2) s3 = simpsonInt(100,1,2) % 6c) g1 = log(2) e1 = ((s1-g1)/g1).*100 e2 = ((s2-g1)/g1).*100 e3 = ((s3-g1)/g1).*100 % 6d) Again i adjusted the value manually. s4 = simpsonInt(16,1,2)
MATLAB/Assignment 2/Newton Method (2)/Newon_script.m
clear format short g x=5; imax=40; tol=1e-6; f ='function1'; df = 'func1_deriv'; xr = Newton(x,f,df,imax,tol); x = linspace(-pi,pi,100)'; y = feval(f,x); yr = feval(f,xr); figure(1);plot(x,y,xr,yr,'o'); title('x root = 4.0885e-06'); clear format short g x=1; imax=40; tol=1e-6; f ='function2'; df = 'func2_deriv'; xr = Newton(x,f,df,imax,tol); x = linspace(-pi,pi,100)'; y = feval(f,x); yr = feval(f,xr); %figure(2);plot(x,y,xr,yr,'o'); %title('x root = 1.1183');
MATLAB/Assignment 2/Newton Method (2)/Newton.m
function root = Newton(x,f,df,imax,tol) % Input: x initial guess at the root % f function to be evaluated % df derivative of f with respect to x % imax iteration maximum % tol convergence criterion % Output: root position of the root % Newton's algorithm % NB: both f and df are user coded functions – you will need to write functions for both!!! % initialise parameters, i.e. % CHECK: verify that derivative is not zero. for i=1:imax xold=x; fold=feval(f,xold); dfold=feval(df,xold); if dfold == 0 error('Problem') end dx = -fold./dfold; x = xold + dx; if abs(dx)<tol; root = x; break end end end
MATLAB/Assignment 2/Bracketing (Multiple Roots) (4)/Newton method bracketing/Newton.m
function root = Newton(x,f,df,imax,tol) % Input: x initial guess at the root % f function to be evaluated % df derivative of f with respect to x % imax iteration maximum % tol convergence criterion % Output: root position of the root % Newton's algorithm % NB: both f and df are user coded functions – you will need to write functions for both!!! % initialise parameters, i.e. % CHECK: verify that derivative is not zero. for i=1:imax xold=x; fold=feval(f,xold); dfold=feval(df,xold); if dfold == 0 error('Problem') end dx = -fold./dfold; x = xold + dx; if abs(dx)<tol; root = x; break end fprintf('Iteration - %i, \t Error - %g \n',i,dx); end end
MATLAB/Assignment 2/Newton method - intersection of two functions (3)/Newton.m
function root = Newton(x,f,df,imax,tol) % Input: x initial guess at the root % f function to be evaluated % df derivative of f with respect to x % imax iteration maximum % tol convergence criterion % Output: root position of the root % Newton's algorithm % NB: both f and df are user coded functions – you will need to write functions for both!!! % initialise parameters, i.e. % CHECK: verify that derivative is not zero. for i=1:imax xold=x; fold=feval(f,xold); dfold=feval(df,xold); if dfold == 0 error('Problem') end dx = -fold./dfold; x = xold + dx; if abs(dx)<tol; root = x; break end fprintf('Iteration - %i, \t Error - %g \n',i,dx); end end
MATLAB/Assignment 1/Simpsons rule(6)/simpsonInt.m
%6a) function [ simpsonInt ] = simpsonInt(n,a,b) function y = evenTot y = 0; for k = 1:n-1 if (mod(k,2) == 0) y = y + feval(@f,(a+(k.*((b-a)/n)))); end end end function g = oddTot g=0; for k = 1:n-1 if (mod(k,2) == 1) g = g + feval(@f,(a+(k.*((b-a)/n)))); end end end simpsonInt = (1/3).*((b-a)/n).*((feval(@f,a))+(feval(@f,b))+(2.*(evenTot))+(4.*(oddTot))); end
MATLAB/Assignment 1/Lab_1 Q 1-4/sin_scale.m
function y = sin_scale (x,sf) % Input: scalar or vector x, sf scale factor % Output: scalar or vector containing corresponding sine values. % y = sin(x); ML intrinsic sine function y = sin(x./sf); % The ‘.’ before an operator implies a vector/array operation end
MATLAB/Assignment 1/Trapezoidal Method(4)/trap1.m
function I = trap1(f,a,b) %Input, f i sintegrand % a and b ar upper and lower limits %Output, s is trapezoidal rule sum if a>b error('b must be greater than a') end h=(b-a); I=h*(feval(f,a)+feval(f,b))/2; end
MATLAB/Assignment 1/Lab_1 Q 1-4/trap1.m
function int = trap1(f,a,b) %Input - f is the integrand input as a string 'f' % - a and b are upper and lower limits of integration %Output - s is the trapezoidal rule sum, equation 3 if a>b error('b must be greater than a') end h=(b-a); int=h*(feval(f,a)+feval(f,b))/2; end
MATLAB/Assignment 1/lab1_trap2 (5)/trap2.m
function I = trap2(f,a,b) %UNTITLED2 Summary of this function goes here % Detailed explanation goes here if a>b error('b must be greater than a') end N=100; I=((b-a)/N)*(0.5)*(feval(f,a)+2*sum(f,i=1..(N-1))+feval(f,b); end
MATLAB/Assignment 1/Simulation on the Nanoscale Lab 1.pdf
Simulation on the Nanoscale ENG-M03: Assignment 1
General information
In the following you will be asked to complete a set of tasks using MATLAB. To help you
accomplish these tasks many questions will include MATLAB code (blue text), this can be
directly copied and pasted into MATLAB. Each question will carry a mark which is indicated
at the questions end. For each question document your answers, be that a numerical or
graphical output in an associated word document, which will be submitted as part of your
continual assessment at the of the semester.
Lab. 1
The aim of this laboratory is to cover the basic set-up used throughout this and following
MATLAB (ML) exercises. First we will review the script-subroutine format of ML by
considering a simple example. Following this we will construct a suite of functions that able
to differentiate and integrate a user defined function numerically using MATLAB.
Using a script file (I)
In this first task we are simply going to plot a function using a script file.
(a) Open a MATLAB script file, save it to your personnel directory as lab1_sin.m. This
script will be used to call intrinsic and user-defined functions in ML.
(b) Type:
clear - clears all variables from the workspace
format short g - formats the outputted numerical data
(c) MATLAB has many ‘in-built’ functions that can be easily used; this task we will plot
the sine function over the interval [-π π].
(d) Define the vector/array/table x as:
N = 100;
x = linspace(-pi,pi,N)’; or as
x = (-pi:2*pi/(N-1):pi)’;
NB: the semi-colon suppresses output to the command window and the apostrophe
transposes the default column vector to a row vector.
(e) Use the following code to plot sin(x) as a function of x:
y = sin(x);
figure(1), plot(x,y)
[1 mark]
Defining a function (II)
A function or subroutine in any programming language is essentially a black-box operator –
i.e. you give the function some value and then it operates on that value and returns the
result of that operation. The sine function in the previous task is an example of this already
built-in to MATLAB. In this task we will define our own MATLAB function, which will scale the
argument given to the sine function, i.e. ⁄ , where a will be user defined.
(a) Open a MATLAB function file
(b) The function file is markedly different to the script file you previously opened and
will look something like:
function [ output_args ] = Untitled( input_args ) %UNTITLED Summary of this function goes here % Detailed explanation goes here:
End
The function starts with the word function and ends with the word end
Untitled – is the default name of the function – the function name MUST be the same as
the name you save it has!!!
Input arguments into the function are written in the curved brackets (separated by
commas) next the function name.
Output arguments must be defined in both the body of the function and in the square
brackets next to the function declaration.
A simple example of a function: function [ y1,y2 ] = multiply_10_and_100( x1,x2 ) % This function multiplies the first input argument by 10 and the second by 100 then return their values in y1 and y2 respectively.
y1 = x1×10; y2 = x2×100;
end
(c) Save the above function as sin_scale.m – make sure you alter the function name to
reflect this change. If you do not do this MATLAB will underline the function name in
yellow.
(d) Modify the function so that you have 1 output and 2 input arguments, your new
function will look like:
function y = sin_scale (x,sf)
% Input: scalar or vector x, sf scale factor
% Output: scalar or vector containing corresponding sine values.
% y = sin(x); ML intrinsic sine function
y = sin(x./sf);
% The ‘.’ before an operator implies a vector/array operation
end
(e) Modify the script previously used in task 1 to plot sin(x) as a function of x, where we
scale x by 0.1.
sf = 0.1
y = sin_scale(x,sf);
figure(2), plot(x,y)
[1 mark]
(f) Using a ‘for loop’ plot the values of sine for the following values sf = [0.1 1 0.5 2]; of
the scale factor. Use the intrinsic ‘zeros’ (type: ‘help zeros’ in the command
window) function to set-aside space for the values of sin(x) for the three values of
scale factor. Use the following code in your script file to help you.
sf=[0.1 0.5 2]; % Scale factor values in an array named sf
N_sf = numel(sf); % Number of elements within sf (ML intrinsic
% function)
sin_array = zeros(N,N_sf);
%Initialise an array of zeros to store data. % Notice the size of the array is predicted by
%the number of values, N and the number of scale factors to be analysed, N_sf.
for ndx = 1: N_sf
sin_array(:,ndx) = sin_function_sf(x,sf(ndx));
end
figure(3),plot(x,sin_array)
title(['A Sine fn; limits = [-pi pi]; Scale factors = ',num2str(sf)])
N.B. you can name your variables as you like, but be careful not to use a ML intrinsic
function. To check this type ‘help variable_name’ if ‘function is not found’ is displayed
you are okay to use the variable name, else change the name.
[2 marks]
Numerical Differentiation
Using the forward difference technique (FDT) we may differentiate a continuous function
over a specified domain. Here, we will consider differentiating the cosine function over the
domain [-pi pi]. Using the FDT the derivative of a function may be expressed as
i
ii
i h
xfxf xf
)()( )('
1
[1]
where i is the index along the x-direction and hi is the step-size between xi+1 and xi; here, we
assume this is a constant.
(a) Define the following input vector x = linspace(-pi,pi,N)’; where N=100; in a script file.
(b) Calculate y = cos(x) – cosine is another built-in MATLAB function.
(c) Define function that inputs both the x and y vectors.
(d) Using equation 1 estimate the cosine derivative.
(e) Plot your answer against the intrinsic ML function ‘-sin’ over the same interval.
[2 marks]
Numerical Integration
In this section we develop some numerical techniques to solve definite integrals of the form
dxxfI b
a
)( [2]
The evaluation of such integrals is usually called quadrature and we will consider two
techniques, namely, the Trapezoidal and Simpson methods.
1. Trapezoidal Method
The area under the curve f(x) from x =a to x = b is approximated by the area beneath a
straight line drawn between the points f(xa,fa) and f(xb,fb). the area below this line is then an
approximation to the integral I, i.e.
))(( 2
1 abffI
ba [3]
This can be evaluated using the following code:
function int=trap1(f,a,b)
%Input - f is the integrand input as a string 'f'
% - a and b are upper and lower limits of integration
%Output - s is the trapezoidal rule sum, equation 3
if a>b
error(‘b must be greater than a’)
end
h=(b-a);
int=h*(feval(f,a)+feval(f,b))/2;
a) Write a function dxxxI 2/
0
)sin(
and evaluate using trap1.
[2 marks]
b) Calculate the analytic value and compare with the numerical answer, i.e. calculate
the percentage error.
[2 marks]
2. Composite Trapezium Rule
Improvement of the accuracy of this technique can be obtained by simply subdividing the
interval of interest into a large number of smaller intervals. For say n regularly spaced
intervals the integral may be approximated by the following expression:
)2( 2
1 1
1
b
n
i
ian fffxI
[4]
Where nabx n
/)( and i
f is the function evaluated at each of the interior points, i.e.
)( ni
xiaxff .
a) Write a ML function to evaluate equation 4.
b) Calculate the integral in the previous section for 100 intervals and compare with the
analytic value, i.e. determine the percentage error.
c) Create a function ⁄ and evaluate it over the interval , with the following number of intervals 5, 50, 100.
[3 marks]
d) Compare each result with the analytic answer i.e. ln(2), i.e. put in a table.
[1 marks]
e) Determine the number of sub-intervals needed to achieve an accuracy to 6 decimal
places.
[2 mark]
3. Simpson’s Rule
An improvement to the composite trapezoidal rule can be achieved by simply replacing the
straight line segments across each interval by parabolic segments i.e. use Simpson’s rule.
From the lecture notes the approximation to the integral across the interval [a b] can be
expressed as follows:
)24( 3
1 1
1
1
1
b
n
onlyeven i
i
n
onlyodd i
ian ffffxI
[5]
a) Write a ML function to evaluate equation 5.
b) Calculate the value of the integral for ⁄ and evaluate it over the interval
, with the following number of intervals 5, 50, 100.
[3 marks]
c) Compare each result with the analytic answer i.e. ln(2), i.e. put in a table.
[1 marks]
d) Determine the number of sub-intervals needed to achieve accuracy to 6 decimal
places.
[2 mark]
MATLAB/Assignment 2/Simulation on the Nanoscale Lab 2.pdf
Simulation on the Nanoscale ENG-M03: Assignment 2
General information
In the following you will be asked to complete a set of tasks using MATLAB. To help you
accomplish these tasks many questions will include MATLAB code (blue text), this can be
directly copied and pasted into MATLAB. Each question will carry a mark which is indicated
at the questions end. For each question document your answers, be that a numerical or
graphical output in an associated word document, which will be submitted as part of your
continual assessment at the of the semester.
Lab. 2
The aim of this lab is to develop a suite of codes that are able to locate the roots of a user
defined function i.e. 0)( xf . You will write scripts to call the bisection algorithm and the
Newton-Raphson algorithm discussed in lectures to find the roots of 4 functions. Finally, you
will use these two methods to locate the multiple roots associated with a fifth function. You
will need to develop a strategy that allows you to (i) approximate the root positions and (ii)
apply both bisection and Newton methods using the approximate root positions to find the
roots to a specified accuracy.
Root Solving
1. Bisection Method
Suppose there is one root within the interval [a b] and that the function f(x) is continuous in
this interval. By defining ax 1
and bx 3
as the left and right ends of the interval and
2/)( 312
xxx as the mid-point. We can determine which half-interval the root resides
by the following if else construction:
if f1*f2<0
% Root is in left half, so
d=(x2+x1)/2;
f3=f2;
x3=x2;
else
% Root is in right half, so
d=(x2+x3)/2;
f1=f2;
x1=x2;
end
i.e. if f1*f2>0 implies both are either positive or negative and the root lies in the right
interval. In either case there is no crossing between x1 and x2. If f1*f2<0 the f(x) has changed
sign between x1 and x2, and thus the root is in the left half. Convergence can then be
determined by the magnitude of the variable d i.e. if told convergence is met, where
the tolerance, tol = 1e-6 say.
Pseudo-code: Bisection Method
function root = bisect(f,a,b,imax,tol)
% Input: f function to be evaluated
% a,b interval limits
% imax iteration maximum
% tol convergence criterion
% Output: root position of the root
% initialise parameters, i.e.
x1=a;
f1=feval(f,x1);
x3=b;
f3=feval(f,x3);
% CHECK: verify that there is indeed a root in the interval [a b].
if f1*f3>0
error('No root in the specified interval')
end
% Begin iteration loop i.e. for ndx =1:imax or while(ndx <imax)
% x2=(x3+x1)/2;
% f2=feval(f,x2); %etc…
% Find which half interval contains root – i.e. if-else construction
% test for convergence i.e. d < tol
% CHECK: Excessive iterations
(a) Using the above pseudo-code, construct a bisection function in MATLAB.
(b) Create the function ( ) ( ) ( ), over the interval [ ]
(c) Create the function ( ) √ ( ), over the interval [ ],
where R = 10;
(d) Create a script function to call the bisection function, in the call to the function
include the parameters:
i. f – a string variable – the value of the file name of the functions defined in
(b) or (c)
ii. a – lower limit of the integration interval
iii. b – upper limit of the integration interval
iv. tol – accuracy to which the root is to be found – 1e-6
v. imax – the maximum number of iterations – 50
(e) For each of the functions (b), (c) plot the function over the desired interval and
indicate the position of the root with a ‘o’ marker, e.g. figure(1);plot(x,y,’-
‘,x_root,y_root,’o’), also plot out the value of the root on the title of the graph, e.g.
title(str2num(x_root)).
[4 marks]
2. Newton’s Method
From the lecture notes the Newton (Newton-Raphson) Method can give a better estimate of
the position of a root ( ) through evaluation of the equation:
( )
( )
Here everything except (the root) is known (assuming that the derivative of the function
can be computed analytically or numerically), and so we can solve for it. Thus, we can apply
this equation iteratively:
( )
( )
( )
( )
( )
( )
Newton’s method converges quadratically towards the root in question (however, bad initial
guesses of the desired root inevitably result in the Newton Method diverging equally as
quickly!).
Pseudo-code: Newton Method
function root = newton(x,f,df,imax,tol)
% Input: x initial guess at the root
% f function to be evaluated
% df derivative of f with respect to x
% imax iteration maximum
% tol convergence criterion
% Output: root position of the root
%
% Newton's algorithm
% xnew = xold + dx
% dx = -fold(xold)/dfold(xold)
% NB: both f and df are user coded functions – you will need to write functions for both!!!
% initialise parameters, i.e.
xold=x;
fold=feval(f,xold); %etc…
% CHECK: verify that derivative is not zero.
count=0;
while norm(f)>tol
if df == 0
error('Problem’)
end
count=count+1;
% calculate equation 4
fprintf('Iteration - %i,\t Error - %g\n',count,norm(fr))
% CHECK: Excessive iterations
end
(a) Using the above pseudo-code, construct a Newton function in MATLAB.
(b) Create the function ( ) (
) and its corresponding derivative
( ) Initial guess
(c) Create the function ( ) ( ) and its corresponding derivative ( )
Initial guess
(d) Create a script function to call the Newton function, in the call to the function
include the parameters:
i. f – a string variable – the value of the file name of the functions defined in
(b) or (c)
ii. – initial guess at root location
iii. tol – accuracy to which the root is to be found – 1e-6
iv. imax – the maximum number of iterations – 20
(e) For each of the functions (b), (c) plot the function over the desired interval and
indicate the position of the root with a ‘o’ marker, e.g. figure(1);plot(x,y,’-
‘,x_root,y_root,’o’), also plot out the value of the root on the title of the graph, e.g.
title(str2num(x_root)).
[4 marks]
3. Newton’s Method – Intersection of two functions
As an example of a more complicated problem utilizing Newton’s method, we will next solve
for the roots of the function:
( ) (
) ( )
which may be written as
(
) ( )
Hence, a root of equation ( ) will correspond to the points of intersection of the two
functions on the left and right hand sides of equation ( ).
(a) Create a script function to call the Newton function, in the call to the function
include the parameters:
a. f – a string variable – the value of the file name of the functions defined in
(b) or (c)
b. – initial guess at root location
c. tol – accuracy to which the root is to be found – 1e-6
d. imax – the maximum number of iterations – 20
(b) Plot these two functions in equation ( ) for the range . It should be
clear from this plot that there are many such intersection points in this interval; we
will only consider the root around . Note that an initial guess greater than 4.5
would result in the wrong root being found. Indicate the position of the root with a
‘o’ marker, e.g. figure(1);plot(x,y,’-‘,x_root,y_root,’o’), also plot out the value of the
root on the title of the graph, e.g. title(str2num(x_root)).
(c) Use your Newton algorithm to determine the root around 4 to a tolerance of 1e-10,
How many iterations are needed to achieve this accuracy? Please comment.
[4 marks]
4. Bracketing (Multiple Roots)
As can be seen above many functions have multiple roots. There many methods available to
determine all the roots of a function. Here we consider the method of bracketing. This
method basically consists of evaluating a given function over a given interval and then
determining and storing the location of the roots within that interval. These stored values
can then act as seeds for a root finding algorithms you have developed.
(a) Develop a bracketing routine that will approximately determines the intersection
points of equation ( ) for the interval [4 10]. Use x = linspace(4,10,100); to
evaluate the function in the first instance.
(b) Use both the Bisection and Newton routines to determine the position for all roots
in this interval to a tolerance of 1e-6.
(c) Plot equation ( ) for the range , Indicating the position of all roots with
a ‘o’ marker
[8 marks]