Open Educational Resources

Ordinary Differential Equations: Problems

  1. The current code given in the explicit Euler method assumes that h is chosen such that \frac{t_{max}-t_0}{h} is a natural number. Modify the code so that if this is not the case, the solution for the dependent variable at t_{max} is available.
  2. Consider the following IVP:

        \[\frac{\mathrm{d}x}{\mathrm{d}t}=xt+4x\]


    with the initial condition x(0)=1.

     

    • Use the built-in DSolve function to find the exact solution and plot x for t\in[0,1].
    • Use the explicit Euler method with h=0.1 to find the values of x_i. Plot the exact solution overlapping the obtained numerical solution.
    • Use the implicit Euler method with h=0.1 to find the values of x_i. Plot the exact solution overlapping the obtained numerical solution.
    • Use Heun’s method with h=0.1 to find the values of x_i. Plot the exact solution overlapping the obtained numerical solution.
    • Use the midpoint method with h=0.1 to find the values of x_i. Plot the exact solution overlapping the obtained numerical solution.
    • Use the classical Runge-Kutta method with h=0.1 to find the values of x_i. Plot the exact solution overlapping the obtained numerical solution.
  3. Consider a mass-spring system with k=m=1 units, g=10 units with the initial conditions x(0)=x'(0)=0.
    • Use the built-in DSolve function to find the exact solution and plot the position x and velocity x' for t\in[0,2\pi].
    • Use the explicit Euler method with h=0.2 to find the values of x_i and x'_i. Plot the exact solutions overlapping the obtained numerical solutions.
    • Use the classical Runge-Kutta method with h=0.2 to find the values of x_i and x'_i. Plot the exact solution overlapping the obtained numerical solution.
  4. Consider the following IVP:

        \[y''(t)+9y=0\]


    with the initial conditions y(0)=1 and y'(0)=0. Use the implicit Euler method to find a numerical solution to the given IVP for t\in[0,5] with h=0.1. Compare with the exact solution.
  5. Consider the following IVP:

        \[y''(t)+2y'(t)+6y=\cos{(2t)}\]


    with the initial conditions y(0)=1 and y'(0)=0. Use the classical 4th-order Runge-Kutta method to find a numerical solution to the given IVP for t\in[0,10] with h=0.1. Compare with the exact solution.
  6. Consider the following nonlinear system of IVPs:

        \[\begin{split}\frac{\mathrm{d}x}{\mathrm{d}t}&=10(y-x)\\\frac{\mathrm{d}y}{\mathrm{d}t}&=x(28-z)-y\\\frac{\mathrm{d}z}{\mathrm{d}t}&=xy-\frac{8}{3}z\\\end{split}\]


    with the initial conditions x(0)=-10, y(0)=-10, and z(0)=30. Use the classical 4th-order Runge-Kutta method to find a numerical solution to the given IVP for t\in[0,20] with h=0.01. Then, using list plots, plot the three curves showing y versus x, z versus x, and y versus z.
    Note: The above system is an example of a Lorenz system of equations which can have chaotic solutions depending on the parameters used and the initial conditions.
  7. Use the finite difference method to find a numerical solution to the following BVP:

        \[4y''(x)-2y-x=0\]


    with the boundary conditions y(0)=0 and y(6)=0. Use n=6 intervals, and compare with the exact solution.
  8. The deflection y as a function of the position x along the length of a simply supported Euler-Bernoulli beam uniformly loaded is given by:

        \[EI\frac{\mathrm{d}^2y}{\mathrm{d}x^2}=\frac{wLx}{2}-\frac{wx^2}{2}\]


    Where E=200,000 MPa is the Young’s modulus, I=30,000 cm^4 is the moment of inertia, w=15 kN/m is the uniformly distributed load, and L=3 m is the length of the beam. The boundary conditions are given by y(0)=y(L)=0. Use the finite difference method to find a numerical solution to the deflection y using n=5 intervals. Compare with the exact solution.
  9. A disgruntled golfer hits a small stone with all his strength. The stone’s acceleration, affected by the drag of the air and by the gravity force, can be modelled using the following pair of differential equations:

    (1)   \[\frac{du}{dt}=-C_d\left(u\sqrt{u^2+v^2}\right)\]


    (2)   \[\frac{dv}{dt}=-C_d\left(v\sqrt{u^2+v^2}\right)-g\]


    where C_d = 0.056 is the drag coefficient, g=9.81 \text{ m/s}^2 in the gravitational acceleration.
    When the stone (assumed to be spherical) leaves the club face, it has an initial speed of U_0 and an initial angle of \theta_0 with respect to the horizontal plane.
    1. Calculate precisely the horizontal distance traveled before the rock first strikes the ground (L), for an initial speed of U_0=265.7 \text{ km/h} and for an initial angle of \theta_0=36^\circ. Use the Euler Explicit Method and select the time step to obtain a solution with a Relative Error Tolerance of 10^{-5}. Describe in detail in your Main Report the procedure used to find L.

    2. Use the same function (or similar) to find accurately and efficiently, the angle corresponding to the maximum horizontal distance that the rock will travel (\theta_{max}). Report it together with the corresponding distance (L_{max}). Describe the procedure used to find \theta_{max}.

    3. Explain why the angle for the maximum distance is not \theta_{max}=45^\circ, You will find the plotting of your results useful to obtain the explanation.

  10. Your company manufactures coil springs of high quality steel and you are asked to program a process controller based on estimated heat treatment times, because actual temperatures cannot be measured. First, the straight steel wires need to be heated until they become ductile, a process called “annealing”. Then they are wound around a cylinder into their coil shape and, finally, they are tempered to give them strength by immersing them in oil (“quenching”) for rapid cooling. You need to estimate the times of annealing, of quenching and total process time.
    In summary, the process is as follows:

    1. Annealing: Steel wires are heated starting from T_{amb} by radiant heaters at T_{rad} that fully surround the wires, while gaining/loosing some heat by free convection to the air at T_{oven}, until T = T_{anneal} is reached.
    2. Forming: When T = T_{anneal} is reached, the wires are removed from the oven and automatically wound around a calibrated cylinder. The forming takes 50 s, during which the wires are cooling by radiation to the surroundings and by free convection to the air at T_{amb}.
    3. Tempering: Finally, the springs are quenched in oil at T_{oil}, an intense free convection process (no radiation), until T = T_{cool} is reached.

    The equation governing the wire temperature T (assumed constant through the cross-section) is

        \[mc_p\frac{dT}{dt}=\dot{Q}_c+\dot{Q}_r\]


    where m is the mass of the wire, A is the surface area of the wire (lateral surface only, neglect ends), the convective heat gain/loss is:

        \[\dot{Q}_c=\begin{cases}-h_{c_{air}}A(T-T_{oven})\quad\text{during annealing}\\-h_{c_{air}}A(T-T_{amb})\quad\text{during forming}\\-h_{c_{coil}}A(T-T_{oil})\quad\text{during quenching}\end{cases}\]


    and the radiative heat gain/loss is:

        \[\dot{Q}_r=\begin{cases}-\epsilon\sigma A(T^4-T_{rad}^4)\quad\text{during annealing}\\-\epsilon\sigma A(T^4-T_{amb}^4)\quad\text{during forming}\end{cases}\]


    The values for the convection heat transfer coefficient h_c are:

        \[ h_{c_{air}}=5\rm W/(m^2K)\]


        \[ h_{c_{oil}}=250\rm W/(m^2K)\]


    The values for temperatures are:

        \[\begin{aligned}T_{amb}&=21^\circ \rm{C} = 294\rm{K}\\ T_{rad}&=850^\circ \rm{C} = 1123\rm{K}\\ T_{amb}&=250^\circ \rm{C} = 523\rm{K}\\ T_{amb}&=650^\circ \rm{C} = 923\rm{K}\\ T_{amb}&=40^\circ \rm{C} = 313\rm{K}\\ T_{amb}&=45^\circ \rm{C} = 318\rm{K}\\ \end{aligned}\]


    The dimensions and properties for the spring wires are: length l= 0.05\rm m, diameter d= 0.003 \rm m, density \rho=8131\rm kg/m^3, c_p=434 \rm J/(kg K), emissivity \epsilon = 0.32, and the Boltzmann constant is \sigma = 5.67\times 10^{-8}\rm W/m^2K^4.

    1. Calculate precisely the time required to complete the annealing (t_{anneal}) and to finish the entire process (t_{total}) for the conditions described above. Use the Euler Explicit Method and select the time step to obtain a solution with a Relative Error of t_{total} below the Tolerance of 10^{-5}. Plot the temperature history for the final curve.
    Important: Use T_{anneal} and the precisely calculated time of the end of the annealing process (t_{anneal}) as Initial Condition and initial time for the forming process. Use the temperature and time at the end of the forming process as Initial Condition and initial time for the tempering process.

    2. To minimize the energy consumption of the process, you need to reduce the temperature of the radiant heaters as low as possible, while ensuring the process still completes in a reasonable time. Use the same function above (or similar) to find accurately and efficiently, the minimum temperature for the radiant heaters T_{{rad}_{min}} that ensures a total time for completion of less than 250 \rm s (t_{total}\le 250). Describe the procedure used to find T_{{rad}_{min}}.

    3. Explain (give physical reason) why does the process fail if T_{rad}= 953 \rm K?
    (Hint: compare intermediate values)

Leave a Reply

Your email address will not be published. Required fields are marked *