Department of Mechanical Engineering
Numerical Techniques Assignments
Student Name: ZHU YU
Student No: X00133408
Date of Submission: 26/11/2018
Supervisor: Ciaran O Sullivan
Declaration of Originality
I hereby certify, that this report, submitted as a part of the BEng (Hons) Mechanical Engineering,
is entirely the work of the author and that any resources used for the completion of the work
done are fully referenced within the text of this report.
Signature: Date:
Q1:
Plot f(x) = 1 / 48 (693x⁶ - 945x⁴ + 315x² - 15) in GeoGebra to find the starting value
for the roots:
The starting value for the roots are -2, -0.7, -0.3, 0.3, 0.7, 2.
Newton-Raphson iteration implementation in EXCEL:
1
Results are shown in the spreed sheet below:
Roots for 1 / 48 (693x⁶ - 945x⁴ + 315x² - 15)=0 are -0.932469, -0.66120939,
-0.23861919, 0.23861919, 0.66120939, 0.932469.
However, if we keep increasing decimal places of these roots, we can get more
accurate answer i.e. more accurate answer for -0.932469 is -0.932469514203152.
So, error for -0.932469 is 1 10 6
error for -0.66120939 is 1 10 8
2
error for-0.23861919 is 1 10 8
error for 0.932469 is 1 10 6
error for 0.66120939 is 1 10 8
error for0.23861919 is 1 10 8
Q2
Plot f(x) = ℯ^x + 1.138 - x³ and g(x) = 0.28561x⁴ + 0.39546x³ - 1.51255x² + 0.156x
+ 0.72 in GeoGebra to find the starting value for the roots:
(i) Starting value for roots of ℯ^x + 1.138 - x³=0 are 2, 4.
Newton-Raphson iteration implementation in EXCEL:
3
Results are shown in the spreed sheet below:
Roots for ℯ^x + 1.138 - x³=0 are 2.108493611(the error is 1 10 9 )
4.4986933(the error is 1 10 7 )
4
(ii) Starting value for roots of 0.28561x⁴ + 0.39546x³ - 1.51255x² + 0.156x +
0.72=0 are -4, -1, 1.
Newton-Raphson iteration implementation in EXCEL:
Results are shown in the spreed sheet below:
5
Roots for 0.28561x⁴ + 0.39546x³ - 1.51255x² + 0.156x + 0.72=0 are
-3.076923077(the error is 1 10 9 )
-0.615384615(the error is 1 10 9 )
1.153846136（the error is 1 10 9 ）
Q3
Starting value for solution of 2-3.8*cos(x)=y^5-2.4*x and x^3+3*y^2=28.5 is (2,2)
Newton-Raphson iteration implementation in EXCEL:
6
Do eight times iteration to improve the initial guess:
Solution for this set of non-linear equations:
(2.731642044,1.644875573) error: 1 10 9
7
Q4
(a)
Change into coupled pair of first order ODEs:
dx
z
dt
dz
x3 x z
dt
(b)
Plot the displacement and velocity of the particle and display the global error bound:
8
Q5
Parameters for my spring-damp system:
M1=180;
M2=50;
K1=15000;
K2=10000;
B=30;
A=20;
w=2;
Results for the displacements and velocities of M1 and M2:
9
Appendix
Code for Q4
Driver:
%NAME: ZHUYU; STUDENT NUMBER:X00133408
% DEdriverPair- Programme to compute solution of a pair of ODEs
% using the 4th order Runge-Kutta method.
% It calls a METHOD function rk4 which in tern will call the
% DETAILS function rhsPairdetails which contains the right side of
the
% system of ODEs written in the vector form dy/dx=f(x,y)
% where y can be a vector of dependent variables.
% variable names:
% indep - stores independent variable
% state - stores dependent variable vector
%
% indepplot - vector to store indep values produced each step
% stateplot - matrix to store dep values produced each step
clear;
help DEdriverPair;
format long;
%% Input of initial values
fred = 0.5;% initialising DUMMY extra parameter
% Initialise independent variable
indep0=input('Enter initial independent variable value:');
indep=indep0; %initialising the independent variable.
% Initialise dependent variable vector
dep0 = input('Enter initial value of first dependent variable:');
dep1 = input('Enter initial value of second dependent variable:');
state=[dep0 dep1];% initialise the state vector
%% Establishing step size and number of steps
final=input('Enter final value of independent variable you want
calculated:');
step=input('Enter independent variable stepsize:');
localerrorestimate= step^5 % local error for rk4
globalerrorestimate = step^4 % global error for rk4
10
nstep=(final-indep0)/step;% calculating the number of steps to be
taken.
%% intialise matrices for storing answers for plot
indepplot=zeros(nstep,1); % produces a matrix of length nstepX1
stateplot=zeros(nstep,2); % produces a matrix of length nstepX2
globalerrorestimatevector = globalerrorestimate*ones(nstep+1,1);
%% main loop
for i=0:nstep
indepplot(i+1,1)=indep; %x
stateplot(i+1,1)= state(1) ;%y
stateplot(i+1,2) = state(2); % z
% call to rk4 function which does work -the sequence of the inputs
to the
% function is the crucial element
state=rk4(indep,state,step,'rhsPairdetails',fred);
indep=indep+step;% incrementing the independent variable by the
stepsize
end
%% Displaying lists of answers and error estimate
disp('indepplot stateplot globalerrorestimatevector')
answers = [indepplot stateplot globalerrorestimatevector]
%% plotting answers
subplot(121)
plot(indepplot,stateplot(:,1),'*')
grid
xlabel('time')
ylabel('displacement')
title('Displacement of the particle' )
text (indep0, dep0, [' global error is' ' '
num2str(globalerrorestimate)] )
subplot(122)
plot(indepplot,stateplot(:,2),'.')
11
grid
xlabel('time')
ylabel('velocity')
title('Velocity of the particle')
Method:
%NAME: ZHUYU; STUDENT NUMBER:X00133408
function Yout=rk4(x,Y, h,derivsRK,param)
% 4th order Runge-Kutta integrator
% Input arguments -
% x=independent variable
% Y=current value of vector of dependent variables
% h=step size
% derivsRK=RHS of the ODE called as derivsRK(x,Y,param)
% param = spare parameter
% Output arguments -
% Yout=new value of Y, vector of dependent variables
%
% Refer to the lecture notes for the meaning of each term
%
hh=0.5*h;
k1=h*feval(derivsRK,x,Y, param);
xhalf=x+hh;
Ytemp= Y+k1/2.;
k2=h*feval(derivsRK,xhalf,Ytemp, param);
Ytemp= Y+k2/2.;
k3=h*feval(derivsRK,xhalf,Ytemp, param);
xfull=x+h;
Ytemp= Y+k3;
k4=h*feval(derivsRK,xfull,Ytemp,param);
Yout= Y+(k1+2.*k2+2.*k3+k4)/6.0;
return;
Details:
%NAME: ZHUYU; STUDENT NUMBER:X00133408
function rhsderiv=rhsPairdetails(w,s, q)
% this function stores the details of the right hand side of a SET
of ODEs
% of the form ds/dw = rhsderiv which is to to be solved numerically,
% where:
% w = independent variable
% s = vector of dependent variables
% q = to allow for the passing of extra parameters
% rhsderiv, contains the details of the right hand side of the
12
SET of ODEs
rhsderiv=[s(2),-s(1)^3-s(1)-s(2)];
return;
Code for Q5
Driver:
%NAME: ZHUYU; STUDENT NUMBER:X00133408
% DEdriverPair- Programme to compute solution of a pair of ODEs
% using the 4th order Runge-Kutta method.
% It calls a METHOD function rk4 which in tern will call the
% DETAILS function rhsPairdetails which contains the right side of
the
% system of ODEs written in the vector form dy/dx=f(x,y)
% where y can be a vector of dependent variables.
% variable names:
% indep - stores independent variable
% state - stores dependent variable vector
%
% indepplot - vector to store indep values produced each step
% stateplot - matrix to store dep values produced each step
clear;
help DEdriverPair;
format long;
%% Input of initial values
fred = 0.5;% initialising DUMMY extra parameter
% Initialise independent variable
indep0=input('Enter initial time:');
indep=indep0; %initialising the independent variable.
% Initialise dependent variable vector
dep0 = input('Enter initial value of x1:');
dep1 = input('Enter initial value of x2:');
dep2 = input('Enter initial value of v1:');
dep3 = input('Enter initial value of v2:');
state=[dep0 dep1 dep2 dep3];% initialise the state vector
13
%% Establishing step size and number of steps
final=input('Enter final value of time you want to calculate:');
step=input('Enter stepsize of time:');
localerrorestimate= step^5 % local error for rk4
globalerrorestimate = step^4 % global error for rk4
nstep=(final-indep0)/step;% calculating the number of steps to be
taken.
%% intialise matrices for storing answers for plot
indepplot=zeros(nstep,1); % produces a matrix of length nstepX1
stateplot=zeros(nstep,4); % produces a matrix of length nstepX2
globalerrorestimatevector = globalerrorestimate*ones(nstep+1,1);
%% main loop
for i=0:nstep
indepplot(i+1,1)=indep; %t
stateplot(i+1,1)= state(1) ;%x1
stateplot(i+1,2) = state(2); %x2
stateplot(i+1,3)= state(3) ;%v1
stateplot(i+1,4) = state(4); %v2
% call to rk4 function which does work -the sequence of the inputs
to the
% function is the crucial element
state=rk4(indep,state,step,'rhsPairdetails',fred);
indep=indep+step;% incrementing the independent variable by the
stepsize
end
%% Displaying lists of answers and error estimate
disp('indepplot stateplot globalerrorestimatevector')
answers = [indepplot stateplot globalerrorestimatevector]
%% plotting answers
subplot(221)
plot(indepplot,stateplot(:,1),'*')
grid
xlabel('t')
ylabel('x1')
14
title('Displacement of M1' )
text (indep0, dep0, [' global error is' ' '
num2str(globalerrorestimate)] )
subplot(222)
plot(indepplot,stateplot(:,2),'*')
grid
xlabel('t')
ylabel('x2')
title('Displacement of M2' )
subplot(223)
plot(indepplot,stateplot(:,3),'.')
grid
xlabel('t')
ylabel('v1')
title('Velocity of M1')
subplot(224)
plot(indepplot,stateplot(:,4),'.')
grid
xlabel('t')
ylabel('v2')
title('Velocity of M2')
Method:
%NAME: ZHUYU; STUDENT NUMBER:X00133408
function Yout=rk4(x,Y, h,derivsRK,param)
% 4th order Runge-Kutta integrator
% Input arguments -
% x=independent variable
% Y=current value of vector of dependent variables
% h=step size
% derivsRK=RHS of the ODE called as derivsRK(x,Y,param)
% param = spare parameter
% Output arguments -
% Yout=new value of Y, vector of dependent variables
%
% Refer to the lecture notes for the meaning of each term
%
hh=0.5*h;
k1=h*feval(derivsRK,x,Y, param);
xhalf=x+hh;
15
Ytemp= Y+k1/2.;
k2=h*feval(derivsRK,xhalf,Ytemp, param);
Ytemp= Y+k2/2.;
k3=h*feval(derivsRK,xhalf,Ytemp, param);
xfull=x+h;
Ytemp= Y+k3;
k4=h*feval(derivsRK,xfull,Ytemp,param);
Yout= Y+(k1+2.*k2+2.*k3+k4)/6.0;
return;
Details:
%NAME: ZHUYU; STUDENT NUMBER:X00133408
function rhsderiv=rhsPairdetails(w,s, q)
% this function stores the details of the right hand side of a SET
of ODEs
% of the form ds/dw = rhsderiv which is to to be solved numerically,
% where:
% w = independent variable
% s = vector of dependent variables
% q = to allow for the passing of extra parameters
% rhsderiv, contains the details of the right hand side of the
SET of ODEs
M1=180;
M2=50;
K1=15000;
K2=10000;
B=30;
A=20;
r=2;
rhsderiv=[s(3) s(4) 1/M1*(-(K1+K2)*s(1)-B*s(3)+K2*s(2)+B*s(4))
1/M2*(K2*s(1)+B*s(3)-K2*s(2)-B*s(4)+A*sin(r*w))];
return;
16
