Homework 12 Solution

12.1 CONCEPT: Characterizing ODEs(3pts)
For each of the following three ODEs …
i. characterize itsorder(e.g.1st-order,2nd-order,etc.),
ii. characterize it as anIVP(Initial Value Problem) orBVP(Boundary ValueProblem),
iii. write the ODE instandardform.

dy⎛d3y ⎞ d2y
(A) xdx⎜dx31⎟dx2, y11,y'22,y''13
⎝ ⎠
⎛d2t⎞
(B) ty tlog⎜dy2⎟, t21,t'12
⎝
1dy2t


⎠

⎛dy⎞

(C) 1
y
dt t⎜dt⎟, y52
⎝ ⎠





Modeling Factory Production (IVP)
Let’s evaluate the productivity of a factory by predicting the number of machines operating (Y) over the course of two months (60 days). On the first day (t0= 0) you observeY0= 20machines are operating.
Then Y will change over time (dY / dt, fortin days) according to the following 3 factors:
• The owner hopes to buy 3 new machines each day (dY/dt =+3),
• The machines have a natural failure rate (dY/dt = – 0.1Y),
• External fluctuating factors affecting machine purchasing and operation (dY/dt =sin[2t/25]).
Combining all these factors leads to thefirst-orderIVP dY
⎛2t⎞
30.1Ysin



with I.C.
Y(0)20

dt ⎝
45⎟

The goal of your analysis to make a plot of number of machines working each dayY(t), and to specifically predict the number of machines working at the end of 2 months (tend= 60days).
Let’s tackle this two ways on the next page:
◦ By-handcomparing the Euler, Modified Euler (Heun’s), and Midpointmethod,
◦ Writing your ownMATLABcode forjustthe Modified Eulermethod.

12.2 APPLICATION: Solving 1st-order IVP, “By-Hand”(14pts) dY
⎛2t⎞
Solve the production model at right using theEULERmethodover the ranget= [0, 20] days with step sizeΔt= 10 days (i.e.calculateYiatt0= 0,t1= 10, andt2=20).
dt30.1Ysin⎜
Y(0)20
45⎟
Show all your work in table form, like we did in class and seen below, showing how you calculate each
ti+1,Yi+1, and “slope” that you use to get each point, with all values accurate tothreedecimals.

i
ti
Yi(known)
Slopei= …
ti+1
Yi+1= …
0
0
Y0= 20
?
?
Y1= …
1
10
Y1= …
?
?
Y2= …
2
20
Y2= …
DONE
DONE
DONE
a) Repeat the work from (a), this time using theMODIFIED EULER(Heun’s) method to calculateY1,Y2.
b) Repeat the work from (a), this time usingMIDPOINTmethod to calculateY1,Y2.
c) Given the exact solution forY(10) = 30.5795, and Y(20) = 34.8952,calculate:
i. theglobalerror (to three decimals) inY2forjustthe Euler method from part(a).
ii. thelocalerror forjustthe single step going fromt1tot2forjustthe Euler method from part(a).


12.3 CODING: Solving the productivity model with yourownModified Euler methodIN MATLAB(7 pts) Go look in the Carmen HW12 folder. You have been given four files:HW12.m,ModEuler.m,f.mand
Y_exact.m, each described below:
f.moutputs the right-hand side of the given factory production IVP in standard form, dY/dt= f(t,Y) :
function dYdt= f(t,Y)
dYdt = 3 – 0.1*Y + sin(2*pi*t/45);
Donotchange anything in this function. Notice that this function …
i. is already “vectorized” (can accept vector inputstandY, and returns a vector outputdYdt),
ii. is a general function oftandY,
iii. the order of the input variables has theindependentvariable (t) first, anddependentvariable
(Y) second. This isreallycrucial later, when we’re using MATLAB’s built-inODE45solver.
ModEuler.mis the functionyoumust complete. It currently only has the function declaration line:
function [t,Y] = ModEuler(t0,tend,Y0,dt)
% Finish the code here to calculate the two output vectors
% t=[t0 t1 ... tend], Y=[Y0 Y1 ... ] given t0, tend, Y0 and dt.
Y_exact.mis also given to provide the exact, analytic solution toY(t) to compare to the ModEuler method later. It’s really long – don’t worry about how I got the analytic solution from your “real math” differential equations course. Just assume it’s correct. ;)
12.3 continued…
HW12.mis the master script. Its job is to solve for the number of operational machinesY(t) with initial numberY0= 20 over the time ranget= [0, 60] days using three different time steps:Δt= 30 days, 3 days, and 0.3 days. It does this by calling your new functionModEuler.mthree times like so:
[t1,Y1]=ModEuler(t0, tend,Y0,30); %Δt= 30d[t2,Y2]=ModEuler(t0, tend,Y0,3); %Δt= 3 d[t3,Y3]=ModEuler(t0, tend,Y0,0.3); %Δt= 0.3d
It also then calculates the exact (analytic) solutions and the error (difference between exact and computational solutions) over the three different time sets:
Y1x =Y_exact(t1); Err1 = abs(Y1x-Y1); Y2x =Y_exact(t2); Err2 = abs(Y2x-Y2); Y3x =Y_exact(t3); Err3 =abs(Y3x-Y3);
Finally, HW12.m makes and saves two plots so you can compare the results for differentΔt:
▪ PlotHW12acompares the three different computational (ti,Yi) with the exactY(t).
▪ PlotHW12bcompares the computational errors inYifor the three choices ofΔt.

Much of the work has already been done for you. But I still need you to …
a) Complete theModEuler.mfunction so that it uses the Modified Euler method to solve theIVPdY/dt= f(t,Y), where f(t,Y) comes from the “slope” functionf.mI gave you.
Make sureModEuler.mworks with general input variables fort0(initialt),tend(finalt),Y0(initial size), anddt(time stepΔt), like its declaration line shows. The function must output vectors for both the time (t) and the operating machines (Y) at each time-step fromt0totendinclusive.
b) Fill in the couple of missing variable names appropriately (currentlyxxx) inHW12.mso that it works again. Also replace the plottitles with more appropriate text for your plot. Other than that, don’t mess much withHW12.m– it’s already designed to work and make the plot for you just fine. If it gives you any problems, it’s because you’re making compatibility errors in theModEuler.mfunction you’re writing.
c) Run the HW12.m script to generate the two plots. Theyshouldshow your Euler-method solutions more closely approaching the exact solution (inred) asΔtgetssmaller.

SUBMIT EVERYTHING FROM HW12.3 ONLINE:
◦ Both your plots (HW12a.pdf,HW12b.pdf) in .pdfform,
◦ Your functionsModEuler.mandHW12.m. (But don’t submit your f.m or Y_exact.m, since you’re not modifying them from what I already gaveyou).
◦ Your answers to thethree questionsbelow in the commentbox:
i. Report the predicted numer of operating machinesYat t = 60 days for each of thethree simulations.
ii. Report theglobal erroratt= 21 days for just the last two simulations (usingΔt= 3 and 0.3days).
iii. How does the change between the two errors in (ii) compare with ouranalytical expectationbased on the error “order” for the Modified Euler method?Hint: yourΔtwent down by afactorof 10, from 3 to 0.3 days – by what factor did your observed error godown??)