Problem Solving in Engineering with Matlab EG-228 PC Lab PART 1 – Introduction to Matlab
Ob jectives
In this lab you will become familiar with the Matlab environment and practice the use of variables, vectors and matrices.
The command window allows you to interact with the program by writing basic in- structions. The workspace allows you to see which variables are currently defined.
Basic operations
The simplest use of Matlab is to perform operations in the command window. Try the following commands and check your understanding of the hierarchy of operations. Check the lecture notes when necesssary.
>> 2*5 + 4/2^2
>> 8*3^2/4
>> 5*(3 – 6/2)/3 + 7*2
>> (4 – 2)/2/4
For introducing a number multiplied by a power of 10 we can use the following short notation
>> 1e4
This is equivalent to 1 × 104.
>> 8.324e8
This is equivalent to 8.324 × 108. >> 2.1e-5
This is equivalent to 2.1 × 10−5.
When an operation is performed in Matlab, the result is assigned to a default variable called ans. Check the workspace every time you write one of the previous expressions in the Matlab command window.
An expression that has been written is remembered by Matlab. Use the up and down arrows (↑ and ↓) to recover previously written expressions.
1
Variables
Define the following variables in the Matlab command window
>> x = 1.549
>> y = -2.001
>> z = 4.578
Check that the variables appear in the workspace. Perform the following operations using the variables
>> a = x + y
>> b = x*y^2 + z
Note that, if a command is finished with a semicolon, the answer is not shown in the screen. Try the following:
>> c = x*y;
>> d = x + y + a;
The value of a variable can be checked by writing its name in the command window. Try the following:
>> a >> d
A variable can be deleted by using the command clear followed by the name of the variable/s that you want to delete. Try the following command:
>> clear c d
Note that the variables c and d have disappeared from the workspace, so they are no longer available to be used. This means, that the following command will produce an error:
>> a + d
The command clear all deletes all the variables in the workspace.
By default, Matlab shows the result on the screen by using four decimal places. You can change the look of the output by using the command format. Try the following commands:
>> format long
>> b
There are other possibilities. Check all of them by typing on the command window
>> help format
2
Functions
Matlab has a huge number of built-in functions that can be used. Try the following functions:
abs: Absolute value.
>> abs(2.4)
>> abs(-1.894)
>> abs(a)
>> abs(b)
sin, cos, tan: cosinus, sinus and tangent (angle in radians).
>> sin(pi/2)
>> cos(pi/4)
>> tan(3*pi/2)
sind, cosd, tand: cosinus, sinus and tangent (angle in degrees).
>> sind(90)
>> cosd(45)
>> tand(270)
sqrt: square root. >> sqrt(2)
There are thousands of Matlab functions. It is important to remember some of the functions that appear in the lecture notes and in the examples, but it is even more impor- tant to develop an ability to find an appropriate function by using the command lookfor or the the Matlab help. See the following two example:
To find the Matlab function required to compute a logarithm in base 2:
Try the following instruction in the command window
>> lookfor logarithm
Open the Matlab help and search for the word logarithm
The Matlab help is usually much better, as it shows examples of how to use a particular function.
Vectors Row vectors
Row vectors can be defined by using square brackets and introducing all the compo- nents of the vector separated by commas or spaces. Try the following commands:
>> v = [2, 8, 4.1, 9, -1.21]
>> w = [-1 0 1.1]
3
It is also possible to define vectors using the command linspace or the syntax: initial:step:final
Check the lecture notes and try the following examples:
>> v1 = linspace(-1,1,10)
>> v2 = linspace(3,5,5)
>> w1 = [0:0.5:2]
>> w2 = [0:6]
>> w3 = [10:-1:6]
Column vectors
Column vectors can be defined by using square brackets and introducing all the com- ponents of the vector separated by semicolons. Try the following commands:
>> s = [4.3; -0.9; 9; -2; 5.81]
>> t = [0; 1; -2; 8]
A column vector can also be defined by transposing a row vector. Try the following command:
>> s1 = transpose(v1)
Operations with vectors
The most common operations with vectors are the addition, subtraction and the scalar product, see the lecture notes. A list of commands is given below, for each one try to anticipate the result that will be output by Matlab and then, check your answer by writing the instruction in the Matlab command prompt. Note that some instructions will produce an error, it is important to understand why!
>> v + v2
>> w-v1
>> w3-w1
>> v*s
>> v*w >> v*w3
Accessing information from vectors
To access the value of a component of a vector, write the name of the vector followed by brackets containing the required component. Try the following commands:
>> v(2)
>> w3(4)
4
It is possible to access to several components at the same time by using the syntax
initial:step:final
Check the lecture notes and try the following commands:
>> v1(2:4)
>> s(5:-2:1)
It is common to require the maximum/minimum value of a vector and where (in which component) is the maximum/minimum located. Try the following commands:
>> h = max(v)
Returns the maximum value of the vector v in a variable named h.
>> [h,p] = max(v)
Returns the maximum value of the vector v in a variable named h and the component of v containing the maximum is given in p. You van check that v(p)=h.
The same applies to the command min.
The number of components of a vector can be computed using the command lenght. Try the following:
>> s = [4.3; -0.9; 9; -2; 5.81]
>> length(s)
Plots
Matlab allows us to represent functions by using two vectors. The first vector must contain the list of x coordinates and the second vector their images. Try the following:
>> x = [0:5:360];
>> y = sind(x);
>> plot(x,y)
The colour of the line can be changed adding an extra argument to the plot command. Try the following:
>> plot(x,y,’k’)
The type of the line can be changed. Try the following:
>> plot(x,y,’r–’)
A symbol can be introduced at each point. Try the following:
>> plot(x,y,’g-s’)
A complete list of all the available options is given in the help. Try the following:
>> help plot
5
After writing help plot, click on the link doc plot to see more examples.
Note that every time we produce a plot, the previous one disappears. To maintain the previous plot and produce a new plot in the same figure, we will use the command hold on. Close the figure and try the following:
>> z = cosd(x);
>> plot(x,y,’k–’)
>> hold on
>> plot(x,z,’r-’)
When different plots are drawn in the same figure, it is necessary to introduce a legend. Try the following:
>> legend(’sinus’,’cosinus’)
Further information can be added to the plot. For instance: A title
>> title(’Trigonometric functions’)
Atextinthexaxis >> xlabel(’x’)
Atextintheyaxis >> ylabel(’y’)
Matrices
A matrix can be defined by writing the coefficients between square brackets. A space or a comma means a change of column and a semicolon a change of row. Try the following instructions:
>> A = [1, 4, 5; 2 6 -1]
>> B = [3 4; 6 9; -1 0; 2 1]
>> C = [3 5 2; -1 0 -1]
>> D = [3 0 1; 0 12 9; 8 7 2]
Matrices can be defined by using existing vectors and/or matrices. Try the following:
>> x = [1, 5, 0]
>> y = [-1, 5, 2]
>> X = [x; y]
>> Y = [x; A; y]
Note that you will receive an error if the dimensions are not consistent. For instance:
>> Z = [x, A]
6
A matrix can also be defined by multiplying a row vector by a column vector. Try the following:
>> x = [1, 5, 0]
>> y = [-1, 5, 2]
>> A = transpose(x)*y
There are some special matrices than can be built without writing all the components one by one. The most useful functions are zeros, ones and eye. See the following examples:
>> A1 = zeros(3,4)
Builds a matrix with 3 rows and 4 columns filled in with zeros
>> A2 = ones(6,2)
Builds a matrix with 6 rows and 2 columns filled in with zeros
>> A3 = eye(5)
Builds an identity matrix of dimension 5. Note that this is equivalent to eye(5,5)
Operations with matrices
The most common operations with matrices are the addition, subtraction and product, see the lecture notes. Try the following instructions. Note that some instructions will produce an error, it is important that you understand why!
>> A + C
>> D – A4
>> A*C
>> A*A1
>> transpose(B)*transpose(A1)
It is possible to solve linear systems of equations by using the backslash operator. Try the following:
>> r = [1; 6; -1]
>> sol = D\r This is the solution of the system of equations with matrix D and right hand side r.
You can check that D*sol=r.
Accessing information from matrices
To access to the value of a component of a matrix, write the name of the matrix followed by brackets that contains the required component. A component is specified by a row and a column separated by a comma. Try the following commands:
>> A(1,3)
>> B(3,2)
7
It is possible to access to several components at the same time by using the syntax initial:step:final for both rows and columns. Check the lecture notes and try the following commands:
>> B(1:2,1)
>> B(1:2:3,2)
>> B(2,1:2)
>> B(:,1)
>> D([1,3],1)
>> D([1,3],:)
>> D([1,3],[1,3])
It is common to require the maximum/minimum value of a matrix and where (in which component) is the maximum/minimum located. Try the following commands:
>> A = [1, 4, 5; 2 6 -1]
>> max(A) Returns a vector with the maximum value of each row of A.
>> [h,p] = max(A) Returns the maximum value of each row of A in a variable named h and the rows of A containing the maximums is given in p.
>> max(max(A)) Returns the maximum value of A.
The same applies to the command min.
The dimension of a matrix can be computed using the command size. Try the fol- lowing:
>> C = [3 5 2; -1 0 -1]
>> size(C)
8
Problem Solving in Engineering with Matlab EG-228 PC Lab PART 1 – Introduction to Programming 1
Ob jectives
In this lab we will practise the creation of Matlab scripts and functions.
For each exercise it is mandatory to create a sketch following the recommendations given in the lecture notes before trying to perform the implementation in Matlab.
Some exercises contain the
Exercises
1. The displacement of a uniformly loaded cantilever beam attached to a wall at x = 0 and free at x = L is given by
y=−wL4 X4−4X3+6X2, 24E I
where w is the load per unit width of the beam, L is the length of the beam, E is the Young modulus, I is the moment of inertia and X = x/L.
a) Write a Matlab function that, given x, L, w, I and E, returns the displacement y.
⋆⋆ Model answer available in Blackboard (Matlab file cantileverUniform.m)
b) Write a Matlab script that defines the load w = 4000N, the moment of inertia I = 863.65×10−8 m4 and the length of the beam L = 2.25m. Then, it computes the displacement at the mid point of the beam for the three materials listed in Table 1 and writes three messages on the screen with the results.
⋆⋆ Model answer available in Blackboard (Matlab file mainCantilever.m)
Table 1: Young modulus for different materials.
2. For a rectangular channel, the Manning’s formula states that
1 A5/3
Q= u S1/2, nP2/3 0
u
where Q is the discharge, Au is the cross sectional area of the channel and Pu is the wetted perimeter, that are expressed in terms of the breath B and the normal depth yu according to
Au =Byu, Pu =B+2yu 9
Material
E
Steel Aluminium Oak
2.0×1011 Pa 6.9×1010 Pa 1.1×1010 Pa
a) Write a Matlab function that, given the breath, the normal depth, the channel slope and the constant n, outputs the discharge.
b) Write a Matlab script that defines the breath B = 50m, the channel slope S0 = 0.002 and the constant n = 0.04. Then, it computes the discharge for two values of the normal depth, namely yu = 5m and yu = 10m, and writes two messages on the screen with the results.
3. Wind tunnel experiments for a particular aerofoil have resulted in the following formulas for the lift and drag coefficients respectively
CL =4.47×10−5α3 +1.15×10−3α2 +6.66×10−2α+1.02×10−1 CD =5.75×10−6α3 +5.09×10−4α2 +1.08×10−4α+1.25×10−2
where the angle α is expressed in degrees.
a) Write a Matlab function that, given the angle, outputs the lift coefficient.
b) Write a Matlab function that, given the angle, outputs the drag coefficient.
c) Write a Matlab script that defines the angle, α = 2.21 degrees. Then, it computes the lift and the drag coefficients and writes two messages on the screen with the results.
4. A cable of length Lc supports a beam of length Lb so that it is horizontal when the weight W is attached at the beam end, see Figure 1. The principles of statics can be used to show that the tension force T in the cable is given by
LbLcW T=22
D Lb−D
where D is the distance of the cable attachment point to the beam pivot.
Figure 1: Cable supporting a beam with an attached weight.
a) Write a Matlab function that, given the distance D, the cable length Lc, the
beam length Lb and the weight W, outputs the tension T.
b) Write a Matlab script that defines the cable length Lc = 5m and the beam length Lb = 3m. Then, it computes the tension at a point with D = Lb/4 for two different loads, namely W = 400N and W = 600N. The program must write two messages on the screen with the results.
10
5. The moment of inertia of a beam with rectangular section is
IR = 1 bh3 (1)
12
where b is the width of the section and h is its height. For a beam with circular section, the moment of inertia is
IC = 1 πd4 (2) 64
where d denotes the diameter.
a) WriteaMatlabfunctionthat,giventhewidthbandtheheighthofarectangular section and the diameter d of a circular section, outputs both the moment of inertia of a beam with rectangular section IR and with circular section IC.
b) Write a Matlab script that defines the width b = 10mm, the height h = 20mm and the diameter d = 15mm. Then, it computes the moment of inertia of a beam with rectangular section IR and with circular section IC. The program must write two messages on the screen with the results.
11
Problem Solving in Engineering with Matlab EG-228 PC Lab PART 1 – Introduction to Programming 2
Ob jectives
In this lab we will practise the control flow (conditionals) and loops.
For each exercise it is mandatory to create a sketch following the recommendations given in the lecture notes before trying to perform the implementation in Matlab.
Exercises
1. The displacement of a punctually loaded cantilever beam attached to a wall at x = 0 and free at x = L is given by
x∈[0,a]
where w is the punctual load, L is the length of the beam, E is the Young modulus, I is the moment of inertia and a is the point where the punctual load is applied.
a) Write a Matlab function that, given x, L, w, I, E and a, returns the displace- ment y.
⋆⋆ Model answer available in Blackboard (Matlab file cantileverPunctual.m)
b) Write a Matlab script that defines the load w = 4000N, the moment of inertia I = 863.65×10−8 m4, the length of the beam L = 2.25m and the Young modulus E = 2.0 × 1011 Pa. Then, it computes the displacement at the mid point of the beam for two different positions of the applied load, namely a = L/3 and a = 2L/3. The program compares the two values and prints on the screen a message with the maximum of both.
⋆⋆ Model answer available in Blackboard (Matlab file mainCantileverPunctual.m)
2. The displacement of a punctually loaded cantilever beam attached to a wall at x = 0 and free at x = L is given in Equation (3) that appears in the previous Worked Example.
a) Create a Matlab script that defines the load w = 4000N, the moment of inertia I = 863.65 × 10−8 m4, the length of the beam L = 2.25m and the Young modulus E = 2.0 × 1011 Pa. Then, it plots the shape of the undeformed and two deformed configurations, with a = L/3 and a = 2L/3.
⋆⋆ Model answer available in Blackboard (Matlab file plotCantileverPunctual.m) 3. Consider the following Taylor series expansion
S (x, n) = (x − 1) − 1 (x − 1)2 + 1 (x − 1)3 − 1 (x − 1)4 + . . . + (−1)n−1 1 (x − 1)n . 234 n
wx2
−
6EI
(3a−x) if
y= (3)
wa2
−6EI(3x−a) if x∈[a,L]
12
a) Write a Matlab function that, given the value of x and a tolerance tol, returns the number of terms k such that the absolute value of the k-th term of the series is less or equal than tol and the value of S(x,k).
⋆⋆ Model answer available in Blackboard (Matlab file taylorTerm.m)
4. For a rectangular channel, the Manning’s formula states that 1 A5/3
Q= u S1/2, nP2/3 0
u
where Q is the discharge, Au is the cross sectional area of the channel and Pu is the wetted perimeter, that are expressed in terms of the breath B and the normal depth yu according to
Au =Byu, Pu =B+2yu
a) Write a Matlab function that, given the breath, the normal depth, the channel
slope and the constant n, outputs the discharge.
b) Write a Matlab script that defines the breath B = 50m, the channel slope S0 = 0.002 and the constant n = 0.04. Then, it plots the discharge for the normal depth varying from yu = 0m to yu = 10m.
5. The displacement of a punctually loaded beam supported at both ends x = 0 and x = L is given by
−wx(L−a)(L2 −x2 −(L−a)2)
6LEI y=
w(L−a) L − 6LEI L−a(x−a)3 +L2 −(L−a)2x−x3
if x∈[0,a] if x∈[a,L]
where w is the punctual load, L is the length of the beam, E is the Young modulus, I is the moment of inertia and a is the point where the punctual load is applied.
a) Write a Matlab function that, given x, L, w, I, E and a, returns the displace- ment y.
b) Create a Matlab script that defines the load w = 4000N, the moment of inertia I = 863.65 × 10−8 m4, the length of the beam L = 2.25m and the Young modulus E = 2.0 × 1011 Pa. Then, it plots the shape of the undeformed beam and two deformed configurations, with a = L/3 and a = 2L/3.
6. For controlling the properties of some materials, it is useful to know the temperature distribution of a rectangular metal plate of dimensions W × L. The temperature depends on the function
2 2 nπx sinh(nπy/L) w(x, y) = π n=1,3,5,… n sin L sinh(nπW/L)
Assuming that W = 4m and L = 3m:
a) Write a Matlab function that, given an integer number k and the point of the plate where the temperature is required (x, y), return the value of w(x, y) where the sum goes up to n ≤ k.
13
b) Write a Matlab script that defines a point x = 1 and y = 1 and plots the temperature at this point as a function of the number of terms in the summation k ∈ [0, 20]. What do you notice as k increases?
7. The static deflection of a square plate clamped on all four sides and subjected to uniform load can be expressed in terms of the following constants Em that are obtained from the truncation of the following infinite set of equations
where
1 α
aiEi +
bimEm = ci, i = 1,3,5,… m=1,3,5,…
ai =
i
bim=8i πm3 1+m2
tanhαi +
i cosh2 αi
4α ci= i − tanh αi
i2 2−1
and αi = iπ/2.
b51
π3i4 cosh2 αi
a) Write a Matlab functions that, given i, returns the coefficient ai.
b) Write a Matlab functions that, given i and m, returns the coefficient bim.
c) Write a Matlab functions that, given i, returns the coefficient ci.
Taking the first four equations, that is i = 1, 3, 5, 7, and truncating the infinite series also to the fourth term, that is m = 1,3,5,7, we obtain the following system of equations in matrix form
a1 +b11
d) Write a Matlab script that, using the functions defined in a) and b), builds and solves the lienar system of equations given in Equation (4). The program must output the solution of the linear system on the screen.
8. Consider two cylinders of two different materials where one cylinder just fits inside the other. The inner radius of the inner cylinder is a, and its outer radius is b. The inner radius of the outer cylinder is also b, and its outer radius is c. The Young’s modulus and the Possion ratio of the inner cylinder are E1 and ν1, respectively, and those of the outer cylinder are E2 and ν2, respectively. The radial stress σrr hoop stress σθθ and radial displacement ur are given by
σrr,i(r) = Ai + Bi r2
σθθ,i(r) = −Ai + Bi (5) r2
b13 a3 + b33
b15 b35
b17 E1 c1
b37 E3 = c3 (4)
b31
b57 E5 c5 b71 b73 b75 a7 + b77 E7 c7
b53
a5 +b55
14
ur,i(r) = −1 + νi Ai + 1 − νi rBi rEi Ei
where i = 1 refers to the inner cylinder and i = 2 to the outer cylinder.
If the outer surface of the outer cylinder is subjected to a compressive radial dis- placement U0 and the inner surface of the inner cylinder has no radial stress, then the following conditions are used to determine Ai and Bi for i = 1, 2.
σrr,1(a) = 0, σrr,1(b) = σrr,2(b)
ur,1(b) = ur,2(b), ur,2(c) = −U0
Introducing these conditions into the Equations (5), the following system of equations in matrix form is obtained:
1+a2 0 0 0 A10 1 b2 −1 −b2 B1= 0 −(1 + ν1) (1 − ν1)b2 (1 + ν2)E1/E2 −(1 − ν2)b2E1/E2 A2 0
0 0 −(1−ν2) (1−ν2)c2 B2 −U0E2c
a) Write a Matlab function that, given a, b, c, E1, ν1, E2, ν2 and U0 builds and solve the linear system and writes on the screen the values of A1, A2, B1 and B2.
b) Write a Matlab function that, given r, Ai, Bi, Ei and νi, returns the radial stress, the hoop stress and the radial displacement.
c) Write a Matlab script that, given a range of values of r and using the functions defined in a) and b), plots in three separate figures the radial stress, the hoop stress and the radial displacement for both the inner and the outer cylinders, as a function of r.
15
Problem Solving in Engineering with Matlab EG-228 PC Lab PART 2 – Interpolation and approximation
Ob jectives
In this lab we will learn useful Matlab functions to solve Engineering problems related to interpolation and approximation.
For each exercise it is mandatory to create a sketch following the recommendations given in the lecture notes before trying to perform the implementation in Matlab.
Matlab functions for interpolation and approximation Function polyfit
Explanation: polyfit builds a polynomial of degree p using the data in the vectors x and y. The vector of coefficients contains the coefficients of the polynomial that has been built:
p(x)=cpxp +cp−1xp−1 +···+c1x+c0
– If the degree p is equal to n − 1, the polynomial interpolates the data. This means that the polynomial passes exactly through the points (x1, y1), (x2, y2), …, (xn,yn).
– If the degree p is less than n − 1, the polynomial approximates the data in a least squares sense.
Usage: c = polyfit(x,y,p)
– Inputs: Two vectors x = [x1,x2,…,xn] and y = [y1,y2,…,yn] and an integer
p.
– Output: One vector of coefficients c = [cp, cp−1, . . . , c0]
Example 1: Write the following instructions in the command window:
>> x = [1 3.25 5.5 7.75 10];
>> y = [2.2094 4.7976 3.4888 6.2615 5.6587];
>> pol = polyfit(x,y,4)
Builds a polynomial of degree 4 that interpolates the data. The polynomial passes through the data because we have 5 points and the polynomial has 5 coefficients (degree 4).
Example 2: Write the following instructions in the command window:
>> x = [1 3.25 5.5 7.75 10];
>> y = [2.2094 4.7976 3.4888 6.2615 5.6587];
>> pol = polyfit(x,y,2)
16
Builds a polynomial of degree 2 that best fits the data. The polynomial does not pass through the data because we have 5 points and the polynomial has only 3 coefficients (degree 2).
The result is pol = [-0.0325 0.7287 1.7858], which means that the poly- nomial is
p(x) = −0.0325×2 + 0.7287x + 1.7858
Function polyval
Explanation: polyval evaluates a polynomial at a given set of points. Usage: Y = polycal(c,X)
– Inputs: The coefficients of a polynomial c = [cp, cp−1, . . . , c0] and a list of points where the polynomial will be evaluated X = [X1, X2, . . . , Xk].
– Output: One vector containing the value of the polynomial at the points in X, namely Y = [Y1,Y2,…,Yk], where Yi = p(Xi), for i = 1,…,k.
Example 1: Write the following instructions in the command window:
>> x = [1 3.25 5.5 7.75 10];
>> y = [2.2094 4.7976 3.4888 6.2615 5.6587];
>> pol = polyfit(x,y,2)
>> f = polyval(pol,1.2)
Evaluates the polynomial approximating the data x and y, with coefficients pol, at x = 1.2.
Example 2: Write the following instructions in the command window:
>> x = [1 3.25 5.5 7.75 10];
>> y = [2.2094 4.7976 3.4888 6.2615 5.6587];
>> pol = polyfit(x,y,4)
>> X = [9.2, 0.1, -1.001]
>> f = polyval(pol,X)
Evaluates the polynomial interpolating the data x and y, with coefficients pol, at x = 9.2, x = 0.1 and x = −1.001.
Function interp1
Explanation: interp1 evaluates a spline at a given set of points. Usage: s = interp1(x,y,X,option). The word could be:
‘linear’: to produce a linear C0 spline. ‘pchip’: to produce a cubic C1 spline. ‘spline’: to produce a cubic C2 spline.
– Inputs: Two vectors x = [x1,x2,…,xn] and y = [y1,y2,…,yn] containing the data, another vector X = [X1, X2, . . . , Xk] containing the points where the spline will be evaluated and an option that will denote the type of spline.
17
– Output: One vector containing the value of the spline at the points in X, namely Y = [Y1,Y2,…,Yk].
Example 1: Write the following instructions in the command window:
>> x = [1 3.25 5.5 7.75 10];
>> y = [2.2094 4.7976 3.4888 6.2615 5.6587];
>> splx = interp1(x,y,[1.2 8.8],’linear ‘)
The result is the value of the linear C0 spline, built from x and y, at the points x = 1.2 and x = 8.8
Example 2: Write the following instructions in the command window:
>> x = [1 3.25 5.5 7.75 10];
>> y = [2.2094 4.7976 3.4888 6.2615 5.6587];
>> splx = interp1(x,y,[1.2 8.8],’pchip ‘)
The result is the value of the cubic C1 spline, built from x and y, at the points x = 1.2 and x = 8.8
Example 3: Write the following instructions in the command window:
>> x = [1 3.25 5.5 7.75 10];
>> y = [2.2094 4.7976 3.4888 6.2615 5.6587]; >> splx = interp1(x,y,[1.2 8.8],’spline ‘)
The result is the value of the cubic C2 spline, built from x and y, at the points x = 1.2 and x = 8.8
18
Exercises
1. The mast of a sailboat is constructed of an experimental aluminium alloy. Tests were performed to define the relationship between stress and strain and the results are:
a) Write a Matlab script that defines two vectors with the stress and strain data ap- pearing in the previous table. Then, it computes and represents on a figure the polynomial that interpolates the data.
⋆⋆ Model answer available in Blackboard (Matlab file mainSailBoatInterpolation.m)
b) Write a Matlab script that defines two vectors with the stress and strain data ap- pearing in the previous table. Then, it computes and represents on a figure the polynomial of degree 2 that approximates the data.
⋆⋆ Model answer available in Blackboard (Matlab file mainSailBoatApproximation.m)
c) Write a Matlab script that defines two vectors with the stress and strain data ap- pearing in the previous table . Then, it computes two polynomials, of degree 3 and 5, that approximate and interpolate the data respectively and it represents the data by using points and the two polynomials using continuous lines with two different colours. Finally, it produces a cubic C1 spline that passes through the data and plots, in a different figure, the data, the polynomial of degree 5 and the spline.
⋆⋆ Model answer available in Blackboard for (c) and (d) (Matlab file mainSailBoat.m)
d) Add the necessary instructions to print on the screen the predicted strain for a stress
of 2 N/cm2 using the two polynomials and the spline produced in a).
⋆⋆ Model answer available in Blackboard for (c) and (d) (Matlab file mainSailBoat.m)
2. The results on the following table show the velocity u of air flowing past a flat surface measured at several distances y away from the surface
a) Write a Matlab script that defines two vectors with the position and velocity data appearing in the previous table. Then, it computes two polynomials, of degree 1 and 5, that approximate and interpolate the data respectively. Finally, the script should represent the data by using points and the two polynomials using continuous lines with two different colours. Add the necessary instructions to add labels, legend, etc, as in the model answer.
b) Design and implement a strategy to estimate the velocity at a distance y = 0.005m. Which of the two polynomials will provide a better approximation of the velocity at y = 0.005m? Why?
c) Design and implement a strategy to estimate the velocity at a distance y = 0.001m. Which of the two polynomials will provide a better approximation of the velocity at y = 0.001m? Why?
Stress [N/cm2]
1.24
3.59
4.97
5.5
6.9
Strain
0.0005
0.0016
0.0032
0.0055
0.0085
y [m]
0
0.002
0.006
0.012
0.018
0.024
u [m/s]
0
0.287
0.899
1.915
3.048
4.299
19
3. After some high accurate measurements, the value of the density and the kinematic viscosity of air at different temperatures is shown in the following table
T [oC]
ρ [Kg/m3]
ν [m2/s]
0
50
100
150
200
250
300
400
500
1.29
1.09 0.946 0.835 0.746 0.675 0.616 0.525 0.457
1.33 × 10−5 1.79 × 10−5 2.30 × 10−5 2.85 × 10−5 3.45 × 10−5 4.08 × 10−5 4.75 × 10−5 6.20 × 10−5 7.77 × 10−5
a) Write a Matlab script that defines three vectors with the temperature (T), density (ρ) and the kinematic viscosity (ν) appearing in the previous table. Then it builds an approximation of the density and the kinematic viscosity for T = 350oC and T = 550oC. Which interpolation or approximation method is more suitable for each case? Why?
4. The effect of temperature on viscosity can be modelled using the so-called Andrade equation
μ = DeB/T (6)
where μ is the dynamic viscosity (measured in Ns/m2), T if the absolute temperature (measured in K) and D and B are parameters. In order to estimate the parameters D and B, the measurements of the following table will be used
a) Write a Matlab script that defines the data of the previous table and implements an appropriate strategy to estimate the parameters D and B of the Andrade equation. Which interpolation or approximation method would you recommend to use? Why?
Hint:
Taking logarithms, we can rewrite Equation (6) as
log(μ) = log(D) + B (7) T
This suggests that after building the two vectors T and μ with the data in the table, it is necessary to build two vectors R and S defined as follows:
– Each component of R is 1 divided by the corresponding component of T (i.e., Ri = 1/Ti).
– Each component of S is the logarithm of the corresponding component of μ (i.e., Si = log(μi).
Then, the problem consists on finding two constants A = log(D) and B such that
S = A + BR 20
T [K]
0.1
5
10
20
30
40
μ [Ns/m2]
1.787
1.519
1.307
1.002
0.7975
0.6529
where data for S and R is available.
Once the constant A is determined, do not forget to obtain D = exp(A).
b) Use the parameters obtained in a) to give an approximation of the dynamic viscosity for a fluid with temperature T = 80 K.
21
Problem Solving in Engineering with Matlab EG-228 PC Lab Root finding
Ob jectives
In this lab we will learn useful Matlab functions to solve Engineering problems related to root finding.
For each exercise it is mandatory to create a sketch following the recommendations given in the lecture notes before trying to perform the implementation in Matlab.
Matlab functions for root finding Function roots
Explanation: roots finds all the zeros of the polynomial given by p(x)=cpxp +cp−1xp−1 +···+c1x+c0
Usage: z = roots(c)
– Inputs: One vector of coefficients c = [cp, cp−1, . . . , c0].
– Output: A vector containing the roots (zeros) of the polynomial with coeffi- cients defined in c.
Example: Write the following instructions in the command window:
>> pol = [-0.0325 0.7287 1.7858];
>> z = roots(pol)
The result is z = [24.6506 -2.2291], which means that the roots of the poly- nomial
p(x) = −0.0325×2 + 0.7287x + 1.7858
are
z1 = 24.6506 z2 = −2.2291
Check that the result is correct by writing in the command window:
>> polyval(pol,z(1))
>> polyval(pol,z(2))
Note: the function polyval was covered in a previous PC lab.
The procedure is illustrated in Figure 2.
22
Figure 2: Quadratic polynomial and the two zeros computed with the function roots. Function fzero
Explanation: fzero finds a zero of a nonlinear function with a given tolerance. This is more generic than roots because it allows you to find zeros of complex functions that are not polynomials.
The method used is not bisection or the Newton method explained in the class but something more advanced. In any case it follows an iterative procedure similar to the bisection or the Newton method, computing better and better approximations to the zero and it stops when the required accuracy is reached.
Usage when the function that you are working with depends upon one variable only: sol = fzero(nameOfFunction,x0,tol)
– Inputs: The name of a function, an initial approximation of the root x0 and (optionally) a tolerance tol. The tolerance, if used, sets the accuracy that you require.
– Output: A zero of the function with the required tolerance. Example 1: Function with one variable.
First you need a Matlab function that contains the function for which you want to compute the zeros. For instance a function named genericFunction.m:
function f = genericFunction(x) %
% Input:
% x: coordinate of a point %
% Output:
% f(x) = sin(x) + atan(x)*log10(x);
f = sin(x) + atan(x)*log10(x);
Write the following instruction in the command window:
>> rootOfF = fzero(@genericFunction, 1, 1e-4)
Note that we need to use @ for Matlab to understand that genericFunction is a Matlab file.
23
The result is rootOfF = 0.099622910166052, which means that if you evaluate the function
f(x) = sin(x) + arctan(x) log10(x)
at x = 0.099622910166052, you should get a zero. In fact, you should get a
value less than the tolerance (10−4 in this example). Check it by writing
>> genericFunction(rootOfF)
Note that the function might have several zeros. It is important to first plot the function and provide an initial approximation that is close to the result that you are looking for. For this example you can obtain the other two zeros by using
>> rootOfF2 = fzero(@genericFunction, 2, 1e-4)
>> rootOfF3 = fzero(@genericFunction, 6, 1e-4)
The procedure is illustrated in Figure 3.
Figure 3: Generic function and the three zeros computed with the function fzero.
Usage when the function that you are working with depends upon several variables:
sol = fzero(@(x) nameOfFunction(x,y),x0,tol)
Note that the difference is that after the @, we need to specify the variable with respect to which you want to find the zero. This is similar to when you differentiate a function with two variables with respect to one only (the other one is taken as constant).
1. Inputs: The name of a function (specifying the variable), an initial approxi- mation of the root x0 and a tolerance tol. The tolerance, if used, sets the accuracy that you require.
2. Output: A zero of the function with the required tolerance. Example 2: Function with two variables.
24
First you need a Matlab function that contains the function for which you want to compute the zeros. For instance a function named genericFunction2.m:
function f = genericFunction2(x, y) %
% Input:
% x, y: coordinates of a point %
% Output:
% f(x,y) = xˆ2*y − yˆ2*x + 2*cos(y) + 3*sin(x);
f = xˆ2*y − yˆ2*x + 2*cos(y) + 3*sin(x);
Write the following instruction in the command window:
>> rootOfF = fzero(@(x) genericFunction2(x,2), 1, 1e-4)
The result is rootOfF = 1.176601062222207, which means that if you evaluate
the function
at x = 1.176601062222207 and y = 2 you should get a zero. In fact, you should
get a value less than the tolerance (10−4 in this example). Check it by writing >> genericFunction2(rootOfF, 2)
The procedure is illustrated in Figures 4 and 5.
Figure 4: Generic function depending on two variables and the plane z = 0. All points on the intersection are roots or zeros of the function.
Note that the function genericFunction2.m depends on two parameters (x, y) and we have computed here the zero with respect to x when y = 2. This is the same as computing the zero of the function
g(x)=2×2 −4x+2cos(2)+3sin(x)
f(x, y) = x2y − y2x + 2 cos(y) + 3 sin(x)
25
Figure 5: The root that we found is the intersection with the plane y = 2 in red.
26
Exercises
1. The principal stresses of a structure can be determined from the roots of the poly- nomial
where
f(σ)=σ3 −C2σ2 −C1σ−C0 C2 = σx+σy+σz
C =τ2+τ2+τ2−σσ−σσ−σσ 1 xy yz zx xy yz zx
C =σσσ+2τττ−στ2−στ2−στ2 0 x y z xy yz zx x yz y zx z xy
σx, σy and σz are the applied normal stresses and τxy, τyz and τzx are the applied shear stresses.
If the roots of the polynomial f(σ) are σ1, σ2 and σ3, the principal shear stresses are
τ12 = (σ1−σ2)/2 τ23 = (σ2−σ3)/2 τ13 = (σ1−σ3)/2
a) Write a Matlab function that, given the values of σx, σy, σz, τxy, τyz and τzx, returns the three coefficients C0, C1 and C2.
⋆⋆ Model answer available in Blackboard (Matlab file coefficientsStress.m)
b) WriteaMatlabscriptthatdefinesσx =100,σy =−60,σz =80,τxy =−40,τyz =50 and τzx = 70. Then, it computes the three coefficients of the polynomial C0, C1 and C2 using the previous function and, by using the function roots, computes the values of σ1, σ2 and σ3. Finally, it outputs on the screen a message with the value of the principal shear stresses.
⋆⋆ Model answer available in Blackboard (Matlab file mainStress.m)
2. The natural frequencies fi of a uniform cantilever beam are related to the roots βi
of the frequency equation
f(β) = cosh(β) cos(β) + 1
whereβi =(2πfi) EI ,misthemassofthebeam,Listhelengthofthebeam,E
2 mL3
is the Young’s modulus and I is the moment of inertia of the cross section.
a) Write a Matlab function that, given the value of β, computes the value of f(β). ⋆⋆ Model answer available in Blackboard (Matlab file freqCantilever.m)
b) Write a Matlab script to plot the function f(β) for β ∈ [0,5].
⋆⋆ Model answer available in Blackboard (Matlab file plotFreqCantilever.m)
c) Consider a beam of steel of length L = 0.9m, with a rectangular cross section 25mm wide and 2.5mm high. Knowing that the mass density of steel is 7 850Kg/m3 and E = 200GPa, determine the lowest two natural frequencies of the beam and outputs a message on the screen with the two values.
⋆⋆ Model answer available in Blackboard (Matlab file mainFreqCantilever.m) 27
3. In the analysis of a nonuniform flow in an open channel of trapezoidal cross section, the ratio of the depth of the fluid to the height of the energy gradient x is determined from
(1 + c0x)2(x2 − x3) = c1 (8) where c0 and c1 depend on the geometry of the channel and the flow rate.
a) Write a Matlab function that, given the value of x, c0 and c1, computes F(x,c0,c1)=(1+c0x)2(x2 −x3)−c1
b) Write a Matlab script that, using the previous function, determines one value of x that satisfies Equation (8) for c0 = 0.4 and c1 = 0.2.
4. Analuminiumcolumn,subjectedtoaneccentricaxialloadP,supportsacompressive stress given by the so-called secant formula
P ec L P σ(P)=A1+r2sec 2r AE
where A is the cross sectional area of the column, e is the eccentricity of the load, c is the half depth of the column, r is the radius of gyration of the cross section, L is the length of the column and E is the Young’s modulus.
a) Write a Matlab function that, given the values of P, A, e, c, r, L and E computes the value of σ.
b) Determine the maximum load P that the column can carry if the maximum stress is not to exceed 120 × 106Pa. Assume the following data:
A = 25800mm2, e = 85mm, c = 170mm, r = 142mm, L = 7100mm and E = 69 × 109Pa.
5. Bernoulli’s equation for fluid flow in an open channel with a small bump is
Q2 Q2
where Q is the volume rate of flow, g is the gravitational acceleration, b is the width of the channel, h0 is the upstream water level, H is the height of the bump and h is the water level about the bump.
a) Write the necessary Matlab functions and/or scripts to obtain the water level about the bump h for a volume rate flow of Q = 1.2m3/s. Assume the following data:
g = 9.81m/s2, b = 1.8m, h0 = 0.6m and H = 0.0075m.
2gb2h2 +h0 = 2gb2h2 +h+H 0
28
Problem Solving in Engineering with Matlab EG-228 PC Lab PART 2 – Numerical integration
Ob jectives
In this lab we will learn useful Matlab functions to solve Engineering problems related to numerical integration.
For each exercise it is mandatory to create a sketch following the recommendations given in the lecture notes before trying to perform the implementation in Matlab.
Matlab functions for numerical integration Function trapz
Explanation: trapz finds the integral of y with respect to x using the trapezoidal rule.
Usage: I = trapz(x,y)
– Inputs: Two vectors x = [x1,x2,…,xn] and y = [y1,y2,…,yn].
– Output: The integral, using the trapezoidal rule, of the data in the vectors x and y.
Example: Write the following instructions in the command window:
>> x = [1.0000 3.2500 5.5000 7.7500 10.0000];
>> y = [2.2094 4.7976 3.4888 6.2615 5.6587];
>> integ = trapz(x,y)
The result can also be visualised using the following instruction >> area(x,y)
29
Exercises
1. Pressure of a wind was measured at various heights on a vertical wall, as shown in Figure 6.
Figure 6: Pressure of a wind on a vertical wall. The pressure centre is defined as
112
hp(h)dh
0
h = 112
p(h)dh 0
(9)
a) Write a Matlab script to compute the pressure centre using the trapezoidal rule and to output a message on the screen with the result.
⋆⋆ Model answer available in Blackboard (Matlab file mainPressureCentre.m) 2. The shape of a cable hanging with no load other than its own weight is a catenary
curve given by
y(x)=10coshx−20, x∈[0,50] (10) 10
where x and y are the horizontal and vertical coordinates measured in meters.
a) It is desired to hang plastic sheeting from the cable to protect the passersby while the bridge is being repainted. Determine how many square meters of sheeting are required. Assume that the bottom edge of the sheeting is located along the x axis at y = 0.
⋆⋆ Model answer available in Blackboard (Matlab files catenary.m and mainCatenary.m)
b) Compute the integral using the trapezoidal rule with 2, 4, 8, 16 and 32 points. Represent the value of the integral as a function of the number of integration points.
⋆⋆ Model answer available in Blackboard (Matlab files catenary.m and mainCatenary2.m) 30
3. The vertical displacement of the surface of an elastic half-space that carries uniform loading over a circular area of radius a at a given point can be shown to be
where
π/2 0
cos2 θ
G(θ) = (r/a)2 − sin2 θ
w(r) = w0
G(θ)dθ, r ≥ a
and w0 is the displacement at r/a.
a) Write a Matlab function that, given the value of θ and r/a, computes the value
of the function G(θ).
b) Write a Matlab script that, using the trapezoidal rule, determines w/w0 for
r/a = 2.
4. The vertical displacement at a certain point of a beam due to a uniform force involves
the evaluation of the following function
1
z2 1 + (2rz)2dz
0
a) Write the necessary Matlab scripts/functions to compute C(r) for a given value of r.
C(r) =
31
Problem Solving in Engineering with Matlab EG-228 PC Lab PART 2 – Ordinary differential equations
Ob jectives
In this lab we will learn useful Matlab functions to solve Engineering problems related to ordinary differential equations.
For each exercise it is mandatory to create a sketch following the recommendations given in the lecture notes before trying to perform the implementation in Matlab.
Matlab functions for ordinary differential equations Function ode45
Explanation: ode45 solves the ODE given by
dZ
in the interval [a, b].
Usage: Z = ode45(@odeRHS, [a, b], Z0)
– Inputs: Name of the function containing the right hand side of the ODE system
written in the form
dZ = F(x, Z), dx
= odeRHS(x, Z) Z(a) = Z0
dx
interval for solving the ODE [a,b] and the initial condition Z0.
Note that we need to use @ for Matlab to understand that odeRHS is a Matlab
file.
– Output: A structure containing the points (field x) where the solution has been approximated and the values of the solution (field y). A structure is a variable than encapsulates several variables inside.
Example 1: Solve the ordinary differential equation dy = y cos(x),
function rhs = odeFunction(x,y) %
% Input
% x: variable of the derivative
% y: solution variable %
with initial condition
Write a Matlab function named odeFunction.m:
dx
y(0) = 1,
32
% Output:
% rhs : right hand side of the ODE %
rhs = y*cos(x);
Write the following instruction in the command window:
>> sol = ode45(@odeFunction, [0, 10], 1);
The solution is stored in the structure sol, see the output produced for this problem
in Figure 7.
Figure 7: Output of a computation using the function ode45.
Check it by writing in the command window:
>> sol.x These are the points where the an approximated solution has been built.
>> sol.y These are the values of the solution.
Note that there is more information in the structure sol that we are not going to use in this course.
Plot the approximated solution:
>> plot(sol.x, sol.y)
The solution is represented in Figure 8
Figure 8: Approximated solution computed using ode45.
33
As you can see, the solution provided by default is not very accurate. To increase the accuracy you can use ode45 with an extra input to request a desired accuracy
Write the following instruction in the command window:
>> sol = ode45(@odeFunction, [0, 10], 1, odeset(‘RelTol ‘, 1e-6) ); where here we are requesting the solution with a relative error less than 10−6. You
can obviously change the requested accuracy.
The output produced for this problem is shown in Figure 9. Note the increase number of points to guarantee a higher accuracy.
Figure 9: Output of a computation using the function ode45. The solution is represented in Figure 10
Figure 10: Approximated solution computed using ode45 requesting higher accuracy.
Example 2: Solve the ordinary differential equation corresponding to the velocity of
a skydiver in free fall (see the example in the lecture notes) dv = g − k v2,
with initial condition
dt m v(0) = 0,
Write a Matlab function named odeFunctionSkydiver.m: 34
function rhs = odeFunctionSkydiver(t ,v) %
% Input
% t: variable of the derivative (time)
% v: solution variable (velocity) %
% Output:
% rhs : right hand side of the ODE %
% Parameters
k = 0.2028; % Air resistance [Kg/m]
m = 80; % Mass [Kg]
g = 9.81; % Gravity [m/sˆ2]
rhs = g − (k/m)*vˆ2;
It is very important to emphasise that the function odeFunctionSkydiver.m has
two inputs but it is actually using only one, x is not used. It is mandatory to
– define functions corresponding to the RHS of an ODE to have two
inputs as shown in the example (even if one is not used). – respect the order of the inputs.
Write the following instruction in the command window:
>> sol = ode45(@odeFunctionSkydiver, [0, 10], 0);
Check it by writing in the command window:
>> sol.x These are the points where the an approximated solution has been built. >> sol.y These are the values of the solution.
Plot the approximated solution:
>> plot(sol.x, sol.y)
The solution is represented in Figure 11
Figure 11: Approximated velocity of a skydiver in free fall computed using ode45.
35
Example 3: Solve the ordinary differential equation corresponding to the position and velocity of a skydiver in free fall (see the example in the lecture notes). The difference with respect to the previous example is that we are interested in position, so the equation is now second order.
with initial condition
d2x = g − k v2, dt2 m
x(0) = 0, v(0) = dx(0) = 0, dt
The ODE is second order so the first thing to do is to write it as a system of 2 ODEs
of order 1. Define a vector Z of 2 components (because the order of the ODE is 2).
The first component is the solution variable x and the second is its first derivative dx.
dt
Z1 x x Z= Z = dx = v
2 dt
For each component of Z, express its derivative in terms of other components of Z.
For the last component we will use the ODE. – First component:
– Second component:
dZ1 =dx=v=Z2 dt dt
dZ2 = d2x = g − k v2 = g − k Z2 dt dt2 m m
Therefore, the system of first ODEs equivalent to the high-order ODE is
with
and
dZ
Z2 F(x,Z)= g−kZ2
m2 0
= F(x, Z) Z(0) = Z0
dx
Z0= 0
Write a Matlab function named odeFunctionSkydiver2.m:
function rhs = odeFunctionSkydiver2(t ,z) %
% Input
% t: variable of the derivative (time)
% z: solution variable (position and velocity , a vector of two components) %
% Output:
% rhs : right hand side of the ODE %
% Parameters
k = 0.2028; % Air resistance [Kg/m]
m = 80; % Mass [Kg]
g = 9.81; % Gravity [m/sˆ2]
rhs = [z(2); g − (k/m)*z(2)ˆ2];
36
Write the following instruction in the command window:
>> sol = ode45(@odeFunctionSkydiver2, [0, 10], [0, 0]);
Note that now sol.y is a matrix that contains in the first row the value of the position and in the second row the value of the velocity, see the output produced for this problem in Figure 12.
Figure 12: Output of a computation using the function ode45 for a system of ODEs.
Check it by writing in the command window:
>> sol.x These are the points where the an approximated solution has been built. >> sol.y(1,:) These are the values of the position.
>> sol.y(2,:) These are the values of the position.
Plot the approximated position and velocity. The results are shown in Figure 13. >> plot(sol.x, sol.y(1,:))
>> plot(sol.x, sol.y(2,:))
Figure 13: Approximated solution for position and velocity.
37
Exercises
1. The transverse displacement y of an Euler beam is governed by the equation d4y =x(1−x), x∈[0,1]
dx4
In addition, the following conditions are satisfied:
w(0) = 0, dy(0) = 0, dy2(0) = −0.2222, dy3(0) = 0.6333 dx dx2 dx3
a) Transform the high-order differential equation to a system of first-order differ-
ential equations.
Model answer:
Define a vector Z of 4 components (because the order of the ODE is 4). The first
component is the solution variable y, the second is its first derivative dy = yx,
the third is its second derivative d2y = y and so on. dx2 xx
Z1 y Z=Z2=yx
Z3 yxx Z4 yxxx
dx
For each component of Z, express its derivative in terms of other components of Z. For the last component we will use the ODE.
First component:
Second component:
Third component:
Fourth component:
dZ1 = dy = yx = Z2 dx dx
dZ2 = d2y = yxx = Z3 dx dx2
dZ3 = d3y =yxxxx =Z4 dx dx3
dZ4 = d4y =yxxxx =x(1−x) dx dx4
Therefore, the system of first ODEs equivalent to the high-order ODE is
with
and
dZ
= F(x, Z) Z(0) = Z0
dx
Z2 F(x, Z) = Z3
Z4 x(1−x)
0 Z0= 0
−0.2222 0.6333
38
b) Write a Matlab function describing the right hand side of the system of first- order differential equations obtained in a).
⋆⋆ Model answer available in Blackboard (Matlab file eulerBeamODErhs.m)
c) Write a Matlab script to plot an approximation of the transverse displacement
of the beam.
⋆⋆ Model answer available in Blackboard (Matlab file mainEulerBeamODE.m)
2. The governing equations of motion, in terms of the slope φ, of a cantilever beam of length L = 1m undergoing large deformations from a follower force P = 8N at the
free end of the beam is
with the initial conditions
d2φ + P sin(φ + π/2) = 0 ds2
φ(0) = 0, dφ(0) = 3.9614. ds
a) Transform the high-order differential equations to a system of first-order differ- ential equations.
b) Write a Matlab function describing the right hand side of the system of first- order differential equations obtained in a).
c) Write a Matlab script to solve the system of ODEs and display the result in a figure.
3. The incompressible flow field in a laminar boundary layer on a flat plate is given by the solution to the ordinary differential equation
where, at z = 0
2d3f+fd2f=0, z∈[0,10] dz3 dz2
f = 0, df = 0, d2f = 0.332043118 dz dz2
a) Transform the high-order differential equation to a system of first-order differ- ential equations.
b) Write a Matlab function describing the right hand side of the system of first- order differential equations obtained in a).
c) Write a Matlab script to solve the system of ODEs and display, in the same figure, the solution f and its first derivative df .
dz
4. The
contact with a cooler fluid is given by the following coupled system of high-order differential equations
d 3 f d 2 f df 2
dz3 +3fdz2 −2 dz +T=0, z∈[0,8] (11)
d2T+3fPrdT=0, z∈[0,8] (12) dz2 dz
equations describing the natural convection along a heated vertical plate in
39
where f is the stream function, T is the temperature and Pr = 0.7 is the Prandtl number. The initial conditions are
f(0)=0, df(0)=0, d2f(0)=0.68 (13) dz dz2
T(0) = 1, dT (0) = −0.50 (14) dz
a) Transform the system of two high-order differential equations to a system of five first-order differential equations.
b) Write a Matlab function describing the right hand side of the system of first- order differential equations obtained in a).
c) Write a Matlab script to solve the system of the stream function and tempera- ture in the interval z ∈ [0, 8]. Display the results in the same figure.
40