Hints on homework #9 ("week 8") in Maple

sm311o_hwk9.mws,wdj,4-6-98

> with(plots):

Section 9.2 Fluid flow in a semi-infinite bay

> psi:=(n,h,x,y)->sinh(Pi*n*x/h)*sin(Pi*n*y/h):
psi(n,h,x,y);
This steam function is the n-th normal mode in an irrotational incompressible fluid in a semi-infinite bay of width h.

> contourplot(psi(1,2,x,y),x=-10..0,y=0..2,axes=BOXED,view=[-10..0,0..2],
contours=[-1000,-100,-30,-10,-5,-2,-1.5,-1,-.5,-.1]);
The contour lines indicate the lines of flow of the fluid (for the first normal mode, with width h=2, in the region -10<x<0).

Problem 2. This is not so easy to do in Maple:

> implicitplot(psi(1,2,x,y) = -1,x=-10..0,y=0..2);
NOTE: First normal mode plots okay.

> implicitplot(psi(2,2,x,y) = -1,x=-10..0,y=0..2);
NOTE: Second normal modes do not plot well using implicitplot or contourplot ( the width does not go to h=2 despite the command y=0..2 ).

To get around this plotting problem, we use the fieldplot command. The velocity vector v of the fluid satisfies v=<M,N>, where M=psi_y, N=-psi_x:

> M:=diff(psi(n,h,x,y),y);
N:=-diff(psi(n,h,x,y),x);
sys:={diff(x(t),t)=subs({x=x(t),y=y(t)},M),diff(y(t),t)=subs({x=x(t),y=y(t)},N)};

> M:=diff(psi(1,2,x,y),y):
N:=-diff(psi(1,2,x,y),x):
fieldplot([M,N],x=-2..0,y=0..2);
Similar to Problem 9, but we have used different values of the parameters.

> M:=diff(psi(2,2,x,y),y):
N:=-diff(psi(2,2,x,y),x):
fieldplot([M,N],x=-1..0,y=0..2);
Problem 2 uses c_2 = -0.001 and h=Pi. We have used c_2=1 and h=2 above.

> M:=diff(psi(3,2,x,y),y):
N:=-diff(psi(3,2,x,y),x):
fieldplot([M,N],x=-1..0,y=0..2);

> M:=diff((0.1)*psi(1,2,x,y)-(0.2)*psi(2,2,x,y),y):
N:=-diff((0.1)*psi(1,2,x,y)-(0.2)*psi(2,2,x,y),x):
fieldplot([M,N],x=-1..0,y=0..2);
This is similar to Problem 5(c), but here h=2 is used instead of h=Pi.

Aside : The above general system of odes for the velocity field of the flow for the n-th normal mode is separable and we can solve them "by hand" or have Maple solve them for y=y(t) as a function of x=x(t). However, it still will not give us a very useful plot because of the complicated nature of the functions involved:

> ode:=diff(y(x),x)=N/M;
Notice that this is a separable but non-linear ode.

> soln:=(n0,h0,x0,y0)->subs({n=n0,h=h0,x=x0},lhs(dsolve(ode,y(x)))):
soln(n,h,x)=c;
y0:=(n,h,x)->solve(soln(n,h,x)=-1,y(x)):
y0(n,h,x);
It's not worth plotting this. The above fieldplots are much better.

Problem 10:

> M:=diff(-psi(1,2,x,y)+(0.2)*x*sin(y*Pi/2),y):
N:=-diff(-psi(1,2,x,y)+(0.2)*x*sin(y*Pi/2),x):
fieldplot([M,N],x=-10..-5,y=0..2,scaling=constrained);

> M:=diff(-psi(1,2,x,y)+(0.01)*x*sin(y*Pi/2),y):
N:=-diff(-psi(1,2,x,y)+(0.01)*x*sin(y*Pi/2),x):
fieldplot([M,N],x=-10..-5,y=0..2,scaling=constrained);

> M:=diff(-psi(1,2,x,y)+2*x*sin(y*Pi/2),y):
N:=-diff(-psi(1,2,x,y)+2*x*sin(y*Pi/2),x):
fieldplot([M,N],x=-3..0,y=0..2,scaling=constrained);

>

Section 10.2 Double integrals

> with(linalg):

Problem 1(a)

> Int(Int(x^2-2*y,y=2..3),x=-1..1);
int(int(x^2-2*y,y=2..3),x=-1..1);
Example 10.2.1. The other examples are done similarly.

Problem 2(a) (the rest of Problem 2 and Problem 3(a), 3(c) are similar.

> Int(Int(x*y,x=-1..1),y=0..1);
int(int(x*y,x=-1..1),y=0..1);

Problem 4(a). The plane passing throught the points (1,0,0), (0,1,0), and (0,0,1) is x+y+z=1,

so z=1-x-y. We integrate this over 0<x<1, 0<y<1-x to get the volume of the tetrahedron:

> Int(Int(1-x-y,y=0..1-x),x=0..1);
int(int(1-x-y,y=0..1-x),x=0..1);

> plot3d(1-x-y,y=0..1-x,x=0..1,axes=normal,color=black);

Problem 6(a)

> g:=(u,v)->u-v; h:=(u,v)->u+v;

> P1a:=plot([g(0,v),h(0,v),v=0..1],color=blue):
P1b:=plot([0,v,v=0..1],color=blue,thickness=2):
P2a:=plot([g(u,0),h(u,0),u=0..1]):
P2b:=plot([u,0,u=0..1],thickness=2):
display([P1a,P1b,P2a,P2b],view=[-2..2,-2..2]);

> jacobian([g(u,v),h(u,v)],[u,v]);
NOTE: Maple does not take the determinant.

The Jacobian of the polar coordinates change-of-variables.

> g:=(u,v)->u*cos(v); h:=(u,v)->u*sin(v);

> P1a:=plot([g(1,v),h(1,v),v=0..1],color=blue):
P1b:=plot([1,v,v=0..1],color=blue,thickness=2):
P2a:=plot([g(u,1),h(u,1),u=0..1]):
P2b:=plot([u,1,u=0..1],thickness=2):
display([P1a,P1b,P2a,P2b],view=[-1..1,-1..1]);

> jac:=jacobian([g(u,v),h(u,v)],[u,v]);
J:=simplify(det(jac));

Problem 7(c). The triangle with vertices (0,0), (1,1), and (-1,1) has boundaries y=x, y=-x, and y=1. We use these to determine our limits of integration:

> Int(Int(y,x=-y..y),y=0..1);
int(int(y,x=-y..y),y=0..1);

> plot3d(1,x=-y..y,y=0..1,axes=framed);

Problem 9(b), Problem 8(a) is similar

> f:=(x,y)->x^2;

> g:=(u,v)->u*cos(v);
h:=(u,v)->u*sin(v);

> Int(Int(f(x,y),x),y);

> Int(Int(f(g(r,theta),h(r,theta)),r=(1/sqrt(2))..1),theta=0..(2*Pi));
int(int(f(g(r,theta),h(r,theta)),r=(1/sqrt(2))..1),theta=0..(2*Pi));

>

>