SDM - Homework 4 - Solved

1.    Consider the spring-mass-damper system in Figure 1. Assume p(t) is zero. The parameters m, k, and c are all positive.

(a)     Develop a free body diagram of the system and write down the equation of motion. Notethat the force from the damper is Fd = −cu˙.

(b)    The behavior of the system depends on the values of the parameters m, k, and c. Develop the general solution to the second order differential equation. What are the three cases that we need to consider?

(c)     Write down the solution for each of the three cases assuming x(0) = 1 and ˙x(0) = 0.

(d)    Rewrite this scalar second order differential equation as a system of first order differentialequations

(e)     Write the system of first order differential equations in matrix vector form u˙ = Au, where u = [x,v]T , with v = x˙.

(f)      Use MATLAB to compute the eigenvalues and eigenvectors of A.

(g)     Using the eigenvalues and eigenvectors, write down the homogeneous solution to thesystem of the first order differential equations.

 

Figure 1: Spring-mass-damper system for Problem 1.

2.    The goal of this problem is to characterize the natural frequencies of the spring mass systemin Figure 2, but now including the effects of gravitation.

(a)     Draw free body diagrams of this system and develop the equations of motion.

(b)    Assume m1 = 2 kg, m2 = 1 kg, k1 = 1000 N/m, and k2 = 2000 N/m and calculate the natural frequencies of the system by assuming the solution is the form of trigonometric functions and then form and solve the resulting eigenproblem.

(c)     Now compute the complete solution by developing the particular solution using the method of undetermined coefficients.

 

Figure 2: Spring mass system for Problem 2.

(d)    Assume initial conditions of x1(0) = 1, x2(0) = 1, ˙x1(0) = x˙2(0) = 0. Plot the homogenous, particular, and complete solutions in one graph and comment on the physical meaning of each.

(e)     Now write the two second order differential equations as a system of first order equationsin matrix form.

(f)      Use the eig command in MATLAB to get the eigenvalues and eigenvectors.

(g)     Write down the complete solution using the homogeneous solution, particular solution,and initial conditions.

3.    Consider the one-dimensional rod in Figure 3. The goal of this problem is to solve the following ODE representing heat transfer along the rod with a fixed temperature at one end and the other end being insulated.

ODE:.

 du

                                                                               BC:          u(0) = 10                (1) = 0

dx

Please compute:

(a)     the homogeneous solution

(b)    the particular solution

(c)     the complete solution

(d)    the integration constants by applying the boundary conditions.

 

                                                               u=10                                                                                              du/dx=0

Figure 3: Heated rod for Problem 3.

4.       Consider the pendulum show in Figure 4. Assume we wish to study the natural motion ofthe system, i.e. the term Mt(t) in the figure is zero. The equation of motion of the system is given by

 

where g is the acceleration due to gravity. This is a nonlinear ordinary differential equation for θ(t). However, for small angles, θ, we can make the approximation sin(θ) ≈ θ.

(a)     Write down the equation of the motion of the pendulum assuming the small angle approximation. We have now linearized the equation. Demonstrate that this equation is now linear.

(b)    Using MATLAB, e.g. to generate a plot, give an estimate of θ for which the small angle approximation ceases to be valid, say 5% error.

(c)     Develop the solution of the differential equation with the small angle approximation.

(d)    Assume we are on Earth, l = 0.5 m, θ(0) = 3o, and θ˙(0) = 0. Write down the solution for the motion of the pendulum under the small angle approximation.

 

Figure 4: Pendulum for Problem 15.

5.       For the following functions compute the Fourier series:

• f1(x) = sin2(x)
 
 
 
(

x+π

• f2(x) = x
if if
− π ≤x<0

0≤x< π
.
(a)     Write down the integral formula to find the coefficients of the Fourier Series for eachfunction; use the period −π ≤ x ≤ π.

(b)    In two separate graphs (one corresponding to each function), plot the actual functionand their Fourier series approximation for N = 5 and N = 100, where N is the number of terms kept in the Fourier series.

6.       Another physical problem, based on PDE models, for which we can use the method of separation of variables, is the two-dimensional Laplace equation. This is a model for steady-state heat transfer in a two-dimensional body.

Consider the rectangular domain x ∈ [0,a] and y ∈ [0,b]. We wish to find the function u(x,y) over the domain. The governing PDE is

∇2u = 0

with boundary conditions:

u(x,0) = 0 u(x,b) = 5 u(0,y) = 0 u(a,y) = 0

(a)     Assume the solution u(x,y) = X(x)Y (y) and develop ODEs for X and Y .

(b)    Using the boundary conditions, determine what must be the form of the functions X(x) and Y (y).

(c)     Determine the coefficients and write down the complete solution for u(x,y).

7.       One estimate that can be used for the condition number of a matrix is the ratio of the largestand smallest eigenvalues. In fact, for real symmetric matrices, that turns out to be exactly the condition number.

Consider the matrix

2

A = 1

0
1

2

1
 
(a)     Perform 3 iterations of the Power method to estimate the largest eigenvalue.

(b)    Perform 3 iterations of the Inverse Power method to estimate the smallest eigenvalue.Make sure to do an LU factorization and use LU for each iteration.

(c)     Compute the ratio of your estimates of the largest and smallest eigenvalues.

(d)    Compare with MATLAB’s values for cond and the values from eig.

8.       Compare the required number of iterations in using QR to find the eigenvalues of F with the number needed when F is preprocessed to Hessenberg form. Consider a reasonable test for convergence.

−631 −798

 F =  900



 −28

96
316

400

−450

14

−48
−156

−195

227

−7

24
144

180

−204

12

−24
−36

−45

51 

−1 

14
9.       The rate of convergence of the power method is influenced by the separation between thelargest eigenvalue and the next largest. Compare the performance of the power method in finding the dominant eigenvalue for the following matrices. Use a reasonable test for convergence.

−631 −798

 F =  900



 −28

96
316

400

−450

14

−48
−156

−195

227

−7

24
144

180

−204

12

−24
−36

−45

51 

−1 

14
−101 −174

 E =  136



 840

2016
51

88

−68

−420

−1008
−12

−20

19 105

252
0

0

0

−32

−84
0 

0 

0 

18

46
10.   We wish to find nontrivial (i.e. non-zero) roots of the following equation:

sin(x) = x3

Use the value of f(x) as the error criterion, i.e. |f(x)| < tol, where the tolerance is given. Write a MATLAB script that does the following:

(a)     Uses a bisection function to find the roots of the nonlinear equation. Use a tolerance of1.0 × 10−6 with a maximum number of iterations of 1000. Use the initial interval of 0.5 to 1.

(b)    Uses a Regula Falsi function to find the roots of the nonlinear equation. Use a tolerance of 1.0 × 10−6 with a maximum number of iterations of 1000. Use the initial interval of 0.5 to 1.

(c)     Uses a Newton’s method function to find the roots of the nonlinear equation. Use atolerance of 1.0 × 10−6 with a maximum number of iterations of 1000.

(d)    Makes a plot of the error vs. the number of iterations that includes each of the differentmethods above. Use the semilogy function to make the plot, i.e. logarithmic scale on the error-axis and linear on the iterations-axis.

11.   Write a MATLAB script that uses your Newton method function to find the root of

f(x) = e−0.5x(4 − x) − 2

using initial guesses of x0 = 2, x0 = 6, and x0 = 8. Explain the results you observe.

12.   A classical problem in engineering is the Catenary problem. This problem is an idealizedmodel of the static equilibrium of a rope or chain suspended from 2 points. The solution to the problem is

where a is a constant to computed and  2). If, additionally, you are given the length of the suspended rope, L, the following equation needs to be solved to compute the constant a:

 

Assuming h = 97 and L = 100, write a MATLAB script that uses your Newton solver to compute a. Plot the resulting Catenary curve.

13.   Consider the nonlinear system below for unknowns x1,x2.

x31 + x2 = 1 x1 + x2 = 0

(a)     Write a MATLAB function that uses the Newton-Raphson method to solve a nonlinearsystem of equations. The function should accept a function handle for a the nonlinear system, a function handle for the Jacobian of the nonlinear system, a stopping tolerance, and maximum number of iterations. The function should return the solution and the error estimate for each iteration. You may use the “\” command to solve the linear system. Use the following as the stopping criterion:

 

γ is the stopping tolerance.

(b)    Write a MATLAB script that uses your Newton-Raphson function to solve the nonlinear

system above, using the initial guess (1,0).

√ 

(c)     Now use an initial guess of (1/ 3,1). Explain the results.

14.   The variation in temperature (in Celsius) of the kinematic viscosity (in m2/s) of SAE 30 oil was found to be as follows:

 

We will develop best fits of this data to the function µ = aexp(−b/T).

(a)     Use a linear regression analysis to determine the parameters a and b.

(b)    Alternatively, we can formulate this as a nonlinear regression problem. Formulate theleast squares nonlinear regression problem. Use Newton’s method to solve the problem. Try both full Hessian and the approximate Hessian neglecting second derivatives of the error function, i.e. the Gauss-Newton method.

(c)     Plot the data and both best fit curves on the same the graph. Discuss.

15.   Recall the pendulum problem. The equation of motion of the system is given by

 

where g is the acceleration due to gravity. We can rewrite second order equations as a system of first order equations and apply solution techniques for initial value problems.

(a)     Write a MATLAB script that uses ode45 to solve this system of differential equations for θ(0) = 5o, θ˙(0) = 0. Assume g = 9.81 m/s2 and l = 1 m. Plot the solution, θ, together with the analytic solution for the small angle approximation, sin(θ) ≈ θ.

(b)    Repeat all steps of the previous part, but now with initial conditions θ(0) = 45o, θ˙(0) = 0.

16.   Consider the differential equation

 

This is a stiff differential equation whose solution has a very sharp initial transient followed by a slower mode. Explicit methods will typically have a much more limited range of stability and, therefore, require much smaller timestep sizes.

(a)     Write a MATLAB script that solves this differential equation using ode45 on the interval t = [0,5], with x(0) = 4. How many time steps did ode45 require?

(b)    Repeat using ode23s, a solver more tailored towards stiff differential equations. How many time steps were required?

17.   The motion of an object, in this case a football, subjected to the forces of gravity and airresistance, can be described by the following system of equations:

 

 where the speed of the object is v = p(x˙)2 + (y˙)2, m = 0.43 kg, g = 9.81 m/s2, ρ = 0.080813 lb/ft3, A is a representative area, and c = .05 for a football (see American Journal of Physics – August 2003 – Volume 71, Issue 8, pp. 791-793). Assume the vertical cross-section of a football is elliptical and use this as the representative area; the length of a football is approximately 11.5 in and its circumference is approximately 21 in.

(a)     Using the equations ˙x = ux and ˙y = uy, rewrite this system of second order equations as a system of first order equations so that, for the vector equation ˙y = f(t,y), y = [x,y,ux,uy]T

(b)    Use ode45 to help answer the following question: on a 60 yard field goal attempt, does the ball clear the crossbar of the goal (10 ft. high) with initial conditions of x(0) = y(0) = 0 and ˙x(0) = y˙(0) = 60.0ft/s? Make plots of x(t) vs. t, y(t) vs. t, and y(t) vs. x(t).

(c)     What if ˙x(0) = 70.0ft/s and ˙y(0) = 58.0ft/s? Make plots of x(t) vs. t, y(t) vs. t, and y(t) vs. x(t).

18.   We’ve previously studied heat conduction in one spatial dimension that included conductionand convection (through Newton’s law of cooling) and saw that, using finite differences, we arrive at a linear system of equations. If we now include radiation, we’ll have instead a nonliner system of equations upon discretization. Consider the following model for one-dimensional heat transfer that includes conduction, convection, and radiation.

 

(a)     By hand, discretize the equations, for an arbitrary number of points N, using a central difference method and write down the resulting equations.

(b)    Explain why the resulting system of equations is nonlinear.

(c)     Instead of solving the nonlinear system of equations, we may instead linearize it. Linearize the radiation term about T¯, a reference temperature that is given, using a Taylor series.

(d)    Again, write down the system of equations now using the linearized form of the radiationterm.

(e)     Show now that we have a linear system of equations.

(f)      Take T∞ = 300K, Ta = 300K, Tb = 2000K, h = 100 W/(m2-K), σ is the Stefan-

Boltzmann constant. Choose an appropriate value for T¯. Discretize both the nonlinear and linearized equations. Use Newton’s method to solve the discretized nonlinear equations. How many points for the discretization do you think you need?

(g)     Repeat the previous part for Tb = 700K. Compare with the previous solutions and discuss.

19.   This problem explores how one has to be careful when choosing which finite difference rules touse when discretizing differential equations. Consider an example of the classical convectiondiffusion transport equation, which describes the transport of a fluid due to convection and diffusion:

au0 − νu00 = 1  0 < x < 1 u(0) = 0 u(1) = 0

Depending on the so-called Peclet number , where ∆x is the spacing between nodes in the finite difference approximation, we can get spurious behavior in the numerical solution if we use a central difference method for the convection (first-order derivative) term. In particular, if Pe > 1 (convection-dominated), then a different scheme should be used.

(a)     Develop the exact solution to the differential equation in terms of a and ν.

(b)    By hand, discretize the differential equation, for an arbitrary number of points N, using a central difference approximation for both terms and write down the resulting equations. When the Pe > 1, this scheme will be “unstable”.

(c)     An alternative scheme is based on “upwinding”. By hand, discretize the differentialequation, for an arbitrary number of points N, using a backward difference rule for the first-order derivative and a central difference rule for the second-order derivative. Write down the resulting equations.

(d)    Write a single Matlab m-file that does the following:

i.         Take a = ν = 1 and N = 21. Compute the numerical solution using both discretizations developed previously. You may use the “\” command in MATLAB to solve the resulting linear systems.

ii.       Plot both numerical approximations and the exact solution on the same plot.

iii.     Take a = 500, ν = 1, and N = 21. Recompute the numerical solutions and plot both approximations and the analytical solution on the same plot.

20.   Runge’s function is given by:

 

Generate 3 sets of points from [-1,1] containing: 5 points, 11 points, 21 points.

(a)     Compute the Lagrange interpolant for each set and plot the interpolating polynomialand the actual function on the same plot.

(b)    Construct fits using linear splines for each set and plot the interpolating function andthe actual function on the same plot. Use the MATLAB function interp1 to construct the spline.

(c)     Construct fits using cubic splines (not-a-knot end conditions) for each set and plot theinterpolating function and the actual function on the same plot. Use the MATLAB function interp1 to construct the spline.

21.   The arclength of the curve y = log(x), 1 < x < 2, is given by

 

Estimate the value of this integral using

(a)     Composite Trapezoidal rule to within an estimated error of 0.1%, using an appropriateerror estimate.

(b)    Composite Simpson’s rule to within an estimated error of 0.1%, using an appropriateerror estimate.

(c)     How many intervals were required for each of the two methods?