# 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),

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))

## P = sin(p y) (e^{a x} + A e^{b x} - 1 -A).

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

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));

**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?

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) (e^{a x} + A e^{b x} - 1 -A).

**data**to get a representative picture of the phase portrait.

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