THE WHALE POPULATION PROBLEM

The model.

x = number of blue whales

y = number of fin whales

G(x,y) = rate of growth for the total whale population (sum of the given growth rates)

>    restart:

>    G := 0.05*x*(1-x/150000)-10^(-8)*x*y + 0.08*y*(1-y/400000)-10^(-8)*x*y;

G := .5e-1*x*(1-1/150000*x)-1/50000000*x*y+.8e-1*y*(1-1/400000*y)

>    plot3d( G, x=0..150000, y=0..400000, axes=BOXED );

[Maple Plot]

The maximum occurs at some value of (x,y), with x > 0 and y > 0.  So far so good.

Now determine x and y that optimize the total growth rate, using multivariable unconstrained optimization techniques from MATH 214 ...

>    dGdx := diff( G, x );

dGdx := .5e-1-.6666666666e-6*x-1/50000000*y

>    dGdy := diff( G, y );

dGdy := -1/50000000*x+.8e-1-.4000000000e-6*y

>    sol := solve( { dGdx=0, dGdy=0 }, {x,y} );

sol := {y = 196544.8172, x = 69103.65549}

Alternatively, we may notice that both partial derivatives are linear in x and y, and so the system to be solved to locate the critical point is a linear system.

We can find its solution with linear algebra methods as follows ...

>    with(linalg): 

Warning, the protected names norm and trace have been redefined and unprotected

>    A := matrix( 2, 2, [2*0.05/150000, 2*10^(-8), 2*10^(-8), 2*0.08/400000 ] );

A := matrix([[.6666666667e-6, 1/50000000], [1/50000000, .4000000000e-6]])

>    b := vector( [ 0.05, 0.08 ] );

b := vector([.5e-1, .8e-1])

>    sol := linsolve( A, b );

sol := vector([69103.65548, 196544.8172])

Excellent, we found the same solution!

Interpretation and discussion: see class notes.

Next - a graphical investigation of the sensitivity of the optimal population levels x and y to the intrinsic growth rate for the blue whales, rB,  set to 0.05 in the model above.

Reassign G, and repeat the solution process from above ...

>    G := rB*x*(1-x/150000)-10^(-8)*x*y + 0.08*y*(1-y/400000)-10^(-8)*x*y;

G := rB*x*(1-1/150000*x)-1/50000000*x*y+.8e-1*y*(1-1/400000*y)

>    dGdx := diff( G, x );

dGdx := rB*(1-1/150000*x)-1/150000*rB*x-1/50000000*y

>    dGdy := diff( G, y );

dGdy := -1/50000000*x+.8e-1-.4000000000e-6*y

>    sol := solve( { dGdx=0, dGdy=0 }, {x,y} );

sol := {x = 12000000.*(250.*rB-1.)/(40000.*rB-3.), y = 7850000000./(40000.*rB-3.)*rB}

We'd like to plot the solutions x and y as functions of rB.  But Maple doesn't recognize x and y as such ...

>    x; y;

x

y

But after using the assign command, it does ...

>    assign(sol);

>    x; y;

12000000.*(250.*rB-1.)/(40000.*rB-3.)

7850000000./(40000.*rB-3.)*rB

Now we can plot x and y as functions of rB ...

>    plot( x, rB=0..0.1, -10000..150000, color=black );

[Maple Plot]

>    plot( y, rB=0..0.1, -10000..400000, color=black );

[Maple Plot]

>    plot( x+y, rB=0..0.1, 0..400000, color=black );

[Maple Plot]

The corresponding quantitative sensitivity analysis goes as follows ...

>    sens_x_rB := diff(x,rB) * rB / x;

>    subs( rB=0.05, sens_x_rB );

sens_x_rB := .8333333333e-7*(3000000000./(40000.*rB-3.)-.4800000000e12*(250.*rB-1.)/(40000.*rB-3.)^2)*rB/(250.*rB-1.)*(40000.*rB-3.)

.8545426808e-1

>    sens_y_rB := diff(y,rB) * rB / y:

>    subs( rB=0.05, sens_y_rB );

-.1502253044e-2

>    sens_xplusy_rB := diff(x+y,rB) * rB / (x+y):

>    subs( rB=0.05, sens_xplusy_rB );

.2111791644e-1

Conclusions for this sensitivity analysis: see class notes.

Need some practice on this stuff?

Then perform a graphical and a quantitative sensitivity analysis to investigate the sensitivity of optimal population levels x and y to the intrinsic growth rate for the fin whale, rF, set to 0.08 in the original model.

Next - a graphical sensitivity analysis of the effect of changes in the maximal sustainable population, KB, set to 150000 in the original model, on the optimal populations x and y ...

NOTE THE USE OF THE UNASSIGN COMMAND BEFORE REDEFINING THE FUNCTION G!!!

>    unassign( 'x', 'y' );

>    G := 0.05*x*(1-x/KB)-10^(-8)*x*y + 0.08*y*(1-y/400000)-10^(-8)*x*y;

G := .5e-1*x*(1-x/KB)-1/50000000*x*y+.8e-1*y*(1-1/400000*y)

>    dGdx := diff( G, x );

dGdx := .5e-1-.10*x/KB-1/50000000*y

>    dGdy := diff( G, y );

dGdy := -1/50000000*x+.8e-1-.4000000000e-6*y

>    sol := solve( { dGdx=0, dGdy=0 }, {x,y} );

sol := {y = 2500000.*(KB-8000000.)/(-100000000.+KB), x = -46000000.*KB/(-100000000.+KB)}

>    assign(sol);

>    plot( x, KB=0..300000, 0..300000, color=black );

[Maple Plot]

>    plot( y, KB=0..300000, 0..400000, color=black );

[Maple Plot]

>    plot( x+y, KB=0..300000, 0..400000, color=black );

[Maple Plot]

Now for the corresponding quantitative analysis ...

>    sens_x_KB := diff(x,KB) * KB / x;

>    subs( KB=150000, sens_x_KB );

sens_x_KB := -.2173913043e-7*(-46000000./(-100000000.+KB)+46000000.*KB/(-100000000.+KB)^2)*(-100000000.+KB)

1.001502253

>    sens_y_KB := diff(y,KB) * KB / y:

>    subs( KB=150000, sens_y_KB );

-.1760602686e-1

>    sens_xplusy_KB := diff(x+y,KB) * KB / (x+y):

>    subs( KB=150000, sens_xplusy_KB );

.2474965983

Conclusions from this sensitivity analysis: see class notes.

You can of course also investigate the sensitivity of the results to changes in KF, set to 400000 in the original model.

This is left as an exercise.

One last sensitivity analysis.  

Suppose alphaB and alphaF are always the same, alpha=alphaB=alphaF.  

In that case, how are the model results affected by changes in alpha?  

Is it ever optimal for one species to become extinct?  

If so, what are the conditions, and which of the whale species should go extinct?

>    unassign( 'x', 'y' );

>    G := 0.05*x*(1-x/150000)-alpha*x*y + 0.08*y*(1-y/400000)-alpha*x*y;

G := .5e-1*x*(1-1/150000*x)-2*alpha*x*y+.8e-1*y*(1-1/400000*y)

>    dGdx := diff( G, x );

dGdx := .5e-1-.6666666666e-6*x-2*alpha*y

>    dGdy := diff( G, y );

dGdy := -2*alpha*x+.8e-1-.4000000000e-6*y

>    sol := solve( { dGdx=0, dGdy=0 }, {x,y} );

sol := {x = .2500000000e15*(-1.+8000000.*alpha)/(-3333333333.+.5000000000e23*alpha^2), y = 200000.*(.6250000000e16*alpha-3333333333.)/(-3333333333.+.5000000000e23*alpha^2)}

>    assign(sol);

>    plot( x, alpha=0..10^(-6), 0..150000, color=black );

[Maple Plot]

>    plot( y, alpha=0..10^(-6), 0..400000, color=black );

[Maple Plot]

Conclusions from this sensitivity analysis: see class notes.