ME G535: Advanced Engineering Mathematics
Assignment-2
Topic: Ordinary and partial differential equations
Maximum marks: 32 (Total weightage: 16 %) Submission deadline: 10 December 2022, 5PM
=============================================================
• Plagiarized codes will not be evaluated. Please strictly follow the code of conduct. Any clarification later
with regard to plagiarism will NOT be entertained.
• Plagiarized codes will be given ZERO marks.
• Each question carries 8 marks.
• Do comment each step in your code. It will act as your notes as you write your code in MATLAB.
• Submit each question as “Q#.m” file, where # is the number of question. Compress all the .m files and
submit the .zip file with name as “FirstName SurName.zip” in Quanta.
=============================================================
1. The motion of a damped spring-mass system is given by the following ordinary differential equation:
d2 x dx
m 2
+c + kx = 0, (1)
dt dt
where x is the displacement from equilibrium position, t is the time, m, c and k are the mass, damping
coefficient and spring stiffness of the system. Let m = 20 kg and k = 20 N/m. Consider three values of
the damping coefficient, viz., c = 5 Ns/m (under-damped), c = 40 Ns/m (critically damped) and c = 200
Ns/m (overdamped). The initial velocity is zero and initial displacement is x = 1 m. Solve the above
equation over the time period 0 ≤ t ≤ 15 s. Use Euler’s, Heun’s, modified Euler’s, RK-2 and RK-4
methods and compare the results obtained for c = 5 Ns/m. Plot the displacement versus time for each of
the damping coefficient on the same plot using RK-4 method.
2. Suppose that the position of a falling object is governed by the differential equation:
d2 x c dx
+ − g = 0, (2)
dt2 m dt
where c = 12.5 kg/s is a first-order drag coefficient, m = 70 kg is the mass, g = 9.81 m/s2 is the
gravitational acceleration. Use the shooting method to solve this equation for position and velocity given
the boundary conditions, x(0) = 0 and x(12) = 500.
3. Take the example of elliptic partial differential equation consisting of Dirichlet boundary conditions as
∂2u ∂2u
+ 2 = 0, (3)
∂x2 ∂y
u(x = 0, y) = 100, (4)
u(x = 12, y) = 100, (5)
u(x, y = 0) = 100, (6)
u(x, 0 < y < 12) = 0, (7)
defined on a square domain x ∈ [0, 12], y ∈ [0, 12]. Consider ∆x 6= ∆y.
(a) Write a MATLAB code to solve the above equation using Gauss-Seidel iterations for solving the
linear system of equations.
1
(b) Consider different discretizations: (i) 20 × 30, (ii) 30 × 30, (iii) 40 × 60, (iv) 80 × 120 divisions in X
and Y directions, respectively. Plot the contours and compare the time taken to find the solution in
the different meshes.
4. Consider a one-dimensional heat conduction equation as
∂T ∂2T
=α 2, (8)
∂t ∂x
defined in a domain 0 ≤ x ≤ L, t ≥ 0 and α = 0.01. The initial condition for the problem is given
as T (x, t = 0) = sin(πx/L), where L = 10. The two boundaries of the domain satisfy T (x = 0, t) =
T (x = L, t) = 0. Consider forward difference in time and central difference in space approximation for
the differential equation.
(a) Write a MATLAB code to find the solution of the problem at each time step, selecting the number
of divisions as ∆x = 0.1, i.e., 100 divisions, and time step size as ∆t = 0.1. Do the calculation till
2000 time steps, i.e., time of 200 s, and plot the contour of temperature variation along the length
of the domain with time, i.e., x vs t with the contour of T .
(b) Now, consider ∆t = 0.4 and do the computations and compare the solution to part (a).
(c) Increase ∆t = 0.51 and run till 950 time steps. What do you observe?
**********
