Worksheet for Project #1

Math 481, Spring 2008

2008-02-13

The original differential equation was for the unknown h(t). From q=Pi*a^2*sqrt(2*g*h) we get h = q^2/(2*Pi^2*a^4*g) therefore we may rewrite the differential equation in terms of q.

restart;

with(DEtools): with(plots):

de_h := Pi * (a + (b-a)/H*h(t))^2 * diff(h(t),t) = -Pi*a^2*sqrt(2*g*h(t)) + Q(t);

Note that the following 4 commands are combined into one execution group for the reason that I explained in class.

Also note that I have terminated the first 3 commands with colons. Change them to semicolons if you want to see their outputs.

subs(h(t)=q(t)^2/(2*Pi^2*a^4*g), de_h):

simplify(%):

simplify(%) assuming a>0, q(t)>0:

de := %;

IC := subs(h(0)=H, a=1, b=7, H=10, g=1, q(0) = Pi*a^2*sqrt(2*g*h(0)));

inflow0 := t -> 0;

DE0 := subs(a=1, b=7, H=10, g=1, Q=inflow0, de);

DEplot(DE0, q(t), t=0..60, q=0..15, [IC], arrows=none, linecolor=black, thickness=1);

inflow1 := t -> 7 + 1.5*sin(t/20);

plot(inflow1(t), t=0..600, 0..10); p1a := %:

DE1 := subs(a=1, b=7, H=10, g=1, Q=inflow1 , de);

DEplot(DE1, q(t), t=0..600, q=0..15, [IC], arrows=none, linecolor=blue, thickness=1, stepsize=1); p1b := %:

display([p1a, p1b], labels=["t", "flow rate"], labeldirections=[horizontal, vertical]);

inflow2 := t -> 7 + 2.5*sin(t/20);

plot(inflow2(t), t=0..600, 0..10); p2a := %:

DE2 := subs(a=1, b=7, H=10, g=1, Q=inflow2, de);

DEplot(DE2, q(t), t=0..600, q=0..15, [IC], arrows=none, linecolor=blue, thickness=1, stepsize=1); p2b := %:

display([p2a,p2b], labels=["t", "flow rate"], labeldirections=[horizontal, vertical]);

inflow3 := t -> piecewise(t<60, 2, t<150, 9, 6);

plot(inflow3(t), t=0..600, 0..10); p3a := %:

DE3 := subs(a=1, b=7, H=10, g=1, Q=inflow3, de);

DEplot(DE3, q(t), t=0..600, q=0..15, [IC], arrows=none, linecolor=blue, thickness=1, stepsize=1); p3b := %:

display([p3a,p3b], labels=["t", "flow rate"], labeldirections=[horizontal, vertical]);