MatLab homework - root solving

profilePhlail
matlab_homework_1__2_answers.rar

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

\\en4a\eg637183$\Personal\MATLAB\Using a script file (I).jpg 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)

\\en4a\eg637183$\Personal\MATLAB\Defining a function (II) (e).jpg 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)

\\en4a\eg637183$\Personal\MATLAB\Defining a function (II) (f).jpg

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

\\en4a\eg637183$\Personal\MATLAB\Numerical Differentiation.jpg 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 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]

MATLAB/Assignment 2/Bracketing (Multiple Roots) (4)/Bisection method bracketing/Bisect roots img.jpg

MATLAB/Assignment 2/Bisecting method (1)/Bisect_fig_1.jpg

MATLAB/Assignment 2/Bisecting method (1)/Bisect_fig_2.jpg

MATLAB/Assignment 1/Images/Defining a function (II) (e).jpg

MATLAB/Assignment 1/Images/Defining a function (II) (f).jpg

MATLAB/Assignment 2/Bracketing (Multiple Roots) (4)/Newton method bracketing/Multiple roots plot.jpg

MATLAB/Assignment 2/Newton method - intersection of two functions (3)/Newton.intesect_fig1.jpg

MATLAB/Assignment 2/Newton Method (2)/Newton_fig1.jpg

MATLAB/Assignment 2/Newton Method (2)/Newton_fig2.jpg

MATLAB/Assignment 1/Images/Numerical Differentiation.jpg

MATLAB/Assignment 1/Images/Using a script file (I).jpg