Reza Malek-Madani

Computer Assignment # 1

In this computer assignment you will be reproducing a fundamental result of ocean circulation, the wind-driven circulation of Henry Stommel (1948). In Stommel's model we encounter a system of nonlinear differential equations of the form

(*) x'=f(x, y), y'= g(x, y),

where x and y are the east-west and north-south coordinates of a typical fluid particle in the Atlantic ocean. The functions f and g take into account the balance of forces due to the prevailing trade winds in the Atlantic as well as the force due to Coriolis.

At this stage of the course we are not ready to derive the equations that I will cite shortly; We are able, however, to solve the initial-value problem that results once f and g are given to us.

The Atlantic ocean of Stommel's model consists of a unit square. The functions f and g have two parameters a and b in them that are affected by the strength of the prevailing winds and the Coriolis parameter. We start our program by defining a constant A by

A = (exp(a)-1)/(1-exp(b))

in terms of which we define the stream function P by

P = sin(p y) (ea x + A eb x - 1 -A).

Finally, the f and g of (*) are defined by

f = dP/dy, g = - dP/dx.

At this stage we have not studied what a stream function is so you are not expected to understand the logic behind the previous three mathematical statements. Suffice it to say that by the end of this semester the mathematical process we just went through to get at f and g will be second nature to you.

Now that f and g are known explicitly, we are ready to plot the phase portrait of (*). The following program does just that for the initial data

((0.5,0.3), (0.5, 0.2), (0.5, 0.1), (0.5, 0.01));

I have arbitrarily chosen a=-10, b=1, and tfinal = 2. Here is the program

A=(Exp[a]-1)/(1-Exp[b]);
P = Sin[Pi y]*(Exp[a x]+A Exp[b x] -1 -A);
f[x_, y_]=D[P, y]; g[x_, y_]=-D[P, x];
eqns={x'[t]==f[x[t],y[t]], y'[t]==g[x[t],y[t]]};
data={{0.5,0.4}, {0.5,0.3}, {0.5, 0.2}, {0.5, 0.1}, {0.5, 0.01}};
tfinal=2; a=-10; b=1;
sol=Table[NDSolve[Flatten[{eqns, 
x[0] == data[[i,1]], 
y[0] == data[[i,2]]}],
{x, y}, {t, 0, tfinal}], 
{i, 1, Length[data]}];
soll1 = Flatten[sol, 1];
output1 = ParametricPlot[Evaluate[{x[t], y[t]} /. soll1], {t, 0, tfinal},
DisplayFunction->Identity];
OutPut = Show[output1,DisplayFunction->$DisplayFunction];

1. Here is the output of the program. Which way is the fluid circulating?

 [Graphics:comp1gr1.gif]

2. Change the parameter values of a and b to -1 and 1, respectively. Plot the phase portrait. How does it differ from the one above?

3. Let's go back to a=-10 and b=1. Change tfinal =1. What is the effect of this parameter on the orbits of the fluid particles?

4. Let b=1 and tfinal = 2. Let a vary from -1 to -10 at increments of -1. Plot the graph of the orbit of (0.5, 0.1) for each a. Combine all orbits on the same screen. What is the effect of varying a on the orbit of this particle? 

5. Replace pi by 2*pi in the definition of P, that is,

P = sin(2p y) (ea x + A eb x - 1 -A).

Make the appropriate changes in the above program and generate the phase portrait of this ocean with tfinal = 2, a = -10, b=1. How does it differ from its counterpart in part 1? You may have to change data to get a representative picture of the phase portrait.

6. Replace pi by 3*pi in the definition of P, that is,

P = sin(3p y) (ea x + A eb x - 1 -A).

Make the appropriate changes in the Mathematica program and generate the phase portrait of this ocean with tfinal = 2, a = -10, b=1. How does it differ from the phase portraits in parts 1 and 5?
go to Top