⎝

₁d_y²t_

⎠

⎛_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 (t₀= 0) you observeY₀= 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 (t_end= 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 1^st-order IVP, “By-Hand”(14pts) _dY

^⎛2t^⎞

1. ⎠

   Solve the production model at right using theEULERmethodover the ranget= [0, 20] days with step sizeΔt= 10 days (i.e.calculateY_iatt₀= 0,t₁= 10, andt₂=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

t_i+1,Y_i+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

2. Repeat the work from (a), this time using theMODIFIED EULER(Heun’s) method to calculateY₁,Y₂.

3. Repeat the work from (a), this time usingMIDPOINTmethod to calculateY₁,Y₂.

4. Given the exact solution forY(10) = 30.5795, and Y(20) = 34.8952,calculate:

   1. theglobalerror (to three decimals) inY₂forjustthe Euler method from part(a).

   2. thelocalerror forjustthe single step going fromt₁tot₂forjustthe 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 …

1. is already “vectorized” (can accept vector inputstandY, and returns a vector outputdYdt),

2. is a general function oftandY,

3. 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. 😉