given in the attachment

profilesethrogen
2016_mech2450_assignment_1.pdf

THE UNIVERSITY OF NEWCASTLE: SCHOOL OF ENGINEERING

Assignment Cover Sheet

All assignments are the responsibility of the student. The School of Engineering cannot accept responsibility for lost or unsubmitted assignments. You are therefore advised to keep a copy.

First Name: -------------------------------------------------

Family Name: -------------------------------------------------

Student Number: -------------------------------------------------

Course Code: *** MECH2450 Assignment #1*** (2016)

Date submitted: -------------------------------------------------

Declaration:

I have read and understand the University of Newcastle’s Policy for the Prevention and Detection of Plagiarism Main Policy Document, which is located at http://www.newcastle.edu.au/policy/academic/general/plagiarism.htm

I declare that, to the best of knowledge and belief, this assignment is my own work, all sources have been properly acknowledged, and the assignment contains no plagiarism. This assignment or any part thereof has not previously been submitted for assessment at this, or any other University.

Student Signature: ------------------------------------------------- Date: ------------------------------

b

MECH2450 Engineering Computations 2 Assignment #1 (2016)

Due Date: 2:00 pm, Monday 15th February 2016. Submit your answers during laboratory directly to Dr Goh Eng Yew.

**You must sign and submit the attached declaration sheet.

**Don’t hesitate to seek help if you need it. Help via email is OK at any time.

Q.1 (From Chapter 1 of the Numerical Analysis part of the Text)

The diagram above shows two corridors of widths a and b (a  b) which meet at right angles. The longest rectangular object (desk or piano etc) of width w (w < a) that can be manoeuvred around the corner has a length l, given by:

  22 /111 mmwambl  (1) provided that the object cannot be moved vertically. The value of m can be found from the condition that in the interval [0, (a/b)1/3 ] m is the only root of the equation:

  02 2222342622  wamwabmmwmwb (2) (a) Write a program using MATLAB to solve Eq. (2) by the Regula Falsi method, and hence find m and l. As a test, for a = 50, b = 70, and w = 40, l = 87.4723. Discuss how the initial points should be chosen.

(b) For thin objects like a ladder, w  0. When w = 0 show that the method will not work if the interval [0, 1] is searched. Explain why this happens. Hint: solve Eq. (2) analytically for w = 0.

(25 marks)

Q.2 a) (from Chapter 2 of the Numerical Analysis part of the Text) Write a program that implements

the LU factorisation algorithm 2.3, with partial pivoting, given in the Text and in the supplementary

section with MATLAB versions of the algorithms. Your code should have separate functions for the

wl b

a

factorization, forward and back-substitution steps. Check your code by analysing the problem solved on

page 41 of the Text:

  

  

    

  

  

  

 

12

2

11

213

122

321

3

2

1

x

x

x

b) (From Chapter 3 of the Numerical Analysis part of the Text) Implement the least squares fitting algorithm 3.2 given in the Lecture Notes and in the supplementary section with MATLAB versions of the algorithms. Use your LU factorisation program (from part (a) of this question) to solve the linear equations and a separate function to evaluate the various expressions. Check your program by analysing the problem given on pages 67-68 of the Text.

c) As an example for an application, you are given the following x-y data representing student attendance at lectures in a particular course:

point x (week) y (no of students at lecture)

1 1 64.0 2 2 68.0 3 3 60.0 4 5 55.0 5 6 62.0 6 9 64.0 7 10 52.0 8 11 57.0

Use your program to find the coefficients for a least squares fit to these data using the expansion: 2 3

0 1 2 3( )f x a a x a x a x    Write your program in MATLAB so that it prints the value of ( )f x at each of the x data points. How well does ( )f x fit the data?

(25 marks)

Question 3 (Solution of ODEs) The SIR model – where “S” stands for susceptible, “I” for infected, and “R” for recovered – is one of the simplest mathematical models for the spread of an infectious disease among a population. The equations are:

)()()(

)()(

tbItItaS dt

dI

tItaS dt

dS



 (1)

and )()()( tItSNtR tot 

where t is time, Ntot is the total number of people in the population, a and b are constants. The constant b is related to the period of infection, p, by p = 1/b. Thus b = 1/3 means that the average period of infection is 3 days. The model also assumes that recovered individuals do not get re-infected. Assume that a single infected person with a disease for which a = 0.0001 and b = 1/3 arrives on an island with a population of 10,000, making the total number of people Ntot = 10,001.

(a) Perform one step (with Δt = 1 day) of the numerical integration of these ODE’s using the fourth order Runge-Kutta method (5. 2) to calculate S (1), I(1) and R(1).

(b) By using the ode45 estimate the number of days required for the disease to ‘run its course’, i.e. for I to be < 0.5. (ode45 – intrinsic MATLAB function to solve systems of ODEs.)

(25 marks)

Question 4 (Solution of PDEs) The torsion (twisting) of a long prismatic bar is governed by the following Poisson equation:

2 2

2 2 2

x y

     

   (1)

where is the torsion function whose value is zero at the edges of the bar. To solve (1) for a bar with a rectangular cross-section, of any lengths a in the y-direction and b in the x-direction, the cross-section is to be divided into an equal number (n) internal nodes in both the x- and y- directions. For node (i,j) with i denoting the x-direction, the finite difference approximations to the two derivatives in (1) are:

2 , 1, , 1,

2 2

2i j i j i j i j

xx h

    

   

and 2

, , 1 , , 1

2 2

2i j i j i j i j

yy h

    

    (2)

where hx is the x-distance between adjacent nodes and hy is the y-distance. (a) Show that the use of (2) leads to the following equation for ,i j :

       

2 2

, 1, 1, , 1 , 1 22 2 1 i j i j i j i j i j

a b a ba b

n        

        

 (3)

where 2

2 2

a a

a b  

 and

2

2 2

b b

a b  

 .

(b) Equation (3) leads to the following form suitable for solution by an iterative technique: ( 1) ( )

, ,

k k

i j i j resid    (4)

where the superscript indicates iteration number. The residual term, resid, is given by

y

a x

x

b

       

2 2

( ) ( ) ( ) ( ) ( ) 21, 1, , 1 , 1 ,2 2 1

k k k k k

i j i j i j i j i j

a b a ba b resid

n         

         

 (5)

Write a MATLAB program to solve (4) and (5) for user inputs of n, lengths a and b, and maximum number of iterations. The program is to output the maximum value of . As a test run your program for n = 50, a = 62cm and b = 75.5cm. Restrict maximum number of iterations to 1000 and set convergence tolerance, . Algorithm for Solving Equation (4) in a rectangular area b x a

In: Number of nodes in both directions, n Lengths a and b Maximum number of iterations, max_iter Convergence tolerance, 

Out: phi(n, n) - the approximate steady state torsion function comment: initialise all elements of phi to zero adash = a^2/(a^2+b^2) bdash = b^2/(a^2+b^2) last=adashbdash(a^2+b^2)/(n+1)^2 do for k = 1, max_iter

maxresid = 0.0 do for i = 2, n+1

do for j = 2, n+1 resid = 0.5(phi(i+1,j) + phi(i-1,j))adash + 0.5(phi(i,j+1) + phi(i,j-1))bdash - phi(i,j) + last

if |resid| > maxresid then maxresid = |resid|

endif phi(i,j) = phi(i,j) + resid

enddo enddo if |maxresid| <  then

exit with solution in phi endif

enddo error: maximum number of iterations exceeded

(25 marks)