UC SCI 111, vM+.
Notebook 5

Off[General::spell1];
Off[Break::nofwd];

Mathematica prerequisites.

Solving differential equations exact: DSolve.

Previously we have already solved a differential equation. This was when we gave three alternative, but equivalent, definitions of the exponential function. The equation was

	d((x(t))/d)x=x(t)

This equation can be solved with Mathematica.  You do this by calling DSolve.

DSolve[D[x[t], t] == x[t], x[t], t]

The arguments you feed to DSolve are the equation, the function you look for and the variable of the function. If we look at the answer there is the weird C[1]. With differential equations need to come initial conditions to obtain uniqueness of a solution. This is what the symbol C[1] refers to. All functions that are a constant times (e^t) will do as a solution. That constant is the C[1]. But you might as well call it something else, C[2] or D[1], if you like. To remove one more possible source of confusion note that (E^x) and Exp[x] are one and the same thing. 

Solving a differential equation is sometimes called integrating it, so let us integrate a differential equation with an initial condition.

DSolve[{D[x[t], t] == x[t], x[0] == 1}, x[t], t]

As you expected the constant C[1] has disappeared.

As examples we integrate a few more differential equations.

DSolve[D[x[t], t] == Sqrt[x[t]], x[t], t]

DSolve[D[x[t], t] == (x[t])^2, x[t], t]

Exercise 3.
Solve the following differential equations

	x'(t) = \(TraditionalForm\`\(\(-1\) - 2\ t\ \(x(t)\)\^3\)\/\(3\ t\^2\ \(x(t)\)\^2\)\)
	
	 x(t) tan t + x'(t) = sin 2t
	
When, i.e. for what values of t, are the solutions defined? Plot some of them.

☺ Some more details on Mathematica.

A mistake people often make is 

	DSolve[D[x[t],t]=x[t],x[t],t]
	
This results in the following error message

	Set::"write" :  
        "Tag \!\(D\) in \!\(\\[PartialD]\_\(t\) \(x[t]\)\) is   \
Protected."

After that some more messages follow, but this first one is the fundamental one. It refers to
	
	D[x[t],t]=x[t]

This is an assignment. And you can't assign the predefined command  D to do something special with the exp function. It's like trying to write 

	2+2=5

Here Mathematica reads

	Plus[2,2]=5

Again, you can't assign the predefined command Plus to do something special with 2 and 2. 
(To be honest, this is a lie. You can force Mathematica to do different things when it encounters the combination of a predefined command and a function you defined. But this trickery is only useful in very special circumstances.) 
So, you have to use the equal sign : "==".

You might have noticed that in the help files of Mathematica the input is often written as

	x'[t] == x[t]

instead of

	D[x[t],t]==x[t]

These are equivalent notations. If you have been trying in vain to find the accent ′ on your keyboard here is how you do it (on a US keyboard). First type x on your keyboard then press Ctrl-6. You have just produced a SuperscriptBox. Now fill the superscript with a special character:
	
	\[Prime]

Move one place to the right with the arrow key to get out of the SuperscriptBox. Then continue your formula. There are loads of special characters. Actually, on MS-Windows these characters will be available system-wide. So you can also include them into MS-Word documents.

Exercise 4. 
a) By solving

	DSolve[D[x[t],t]==1/x[t],x[t],t]

by hand explain why you find two expressions for the solution. 
b) Write down a differential equation that will give five expressions for the solution.

Solving differential equations numerically:  NDSolve.

Not all differential equations can be solved in the sense that an explicit expression for the solution can be derived. In fact, most cases are hopeless. Here is just one example that results in horrendous output. It is quite easy to produce a lot of weird messages by using DSolve. Here is one that is particularly puzzling.

DSolve[{D[x[t], t] == -x[t]^2 + x[t] + x[t]^5}, x[t], t]

The same happens with a differential equation that will come up later

(Subscript[g, out])[xx_] := c*xx/(k + xx);
DSolve[{D[x[t], t] == v - (Subscript[g, out])[x[t]]}, x[t], t]

Here we see that we might not get all solutions and what is worse is that ProductLog is a function that refers to, in some way which we will not explain, solutions to another differential equation, that can neither be solved. So, we haven't solved anything.

For these hopeless cases we  have another Mathematica command, that can help us out. There are algorithms that solve these differential equations numerically. The results are not nice expressions such as Exp[t], instead we get a so called InterpolatingFunction. This yields a numerical approximation to the true solution.

The input to NDSolve consists of three parts.
1.	The differential equation (I1) and its initial conditions (I2).
	If you do not specify enough initial conditions you will get an error message!
  	NDSolve can't deal with arbitrary constants, like C[1]. 
2.	The function you look for (I3).
3.	The domain in which the approximate solutions should be valid (I4).
	This domain is entered as in the Plot command.
This recipe is applied in the next example.

sls=NDSolve[
  {x''[t]==-Sin[x[t]],               
x[0] == 0,(x^′)[0] == 1},           
       x[t],                                           
 {t, -Pi, Pi}                                
]

Plot[x[t] /. sls, {t, -Pi, Pi}];

We can also get some idea of ProductLog. It satisfies a differential equation that can be solved by NDSolve, but not by DSolve.

sls = 
      NDSolve[{x'[t] == x[t]/(t (1 + x[t])), x[1] == 1}, 
        x[t], {t,  .0001, 1}]

Plot[x[t] /. sls, {t,  .0001, 1}];

Exercise 5.
Remember that (Subscript[g, out]) is given by

(Subscript[g, out])[xx_] := c*xx/(k + xx);

a) Use NDSolve to solve numerically

	x'(t) = v - (Subscript[g, out])(x(t))
	
Start with v=1, c=1, k=1 and x(0)=0.
b) Experiment with different initial conditions.

Projects.

Continuously stirred tank reactor.

Consider a continuously stirred tank reactor (CSTR) of volume V  in which per unit of time a volume  F  streams in and out (such that the total volume remains constant). Let x  (in mol/unit of volume) be the concentration of a chemical X . Assume that instreaming liquid does not contain X  but that in some other way I  mol of X  is added per unit of time.

Exercise 6.
a) Show that

	(d/dt)(V x) = -F x + I
	
and hence

	(dx/dt) = - (F/V) x + (I/V)

where F / V  is the turnover rate (or transfer coefficient; it is the fraction of the volume that is replaced per unit of time; in a pharmacological context this is often called clearance and defined as the fraction of body liquid that is completely cleared of pharmacon per unit of time).
b)   Determine the steady state and its stability. Express the relaxation time in terms of the parameters.
c) For those who like to do calculations: find an explicit expression for the solution of the initial value problem. 
Hint: consider the deviation from the steady state as a new variable or apply Exercise 2.

Lakes.

Exercise 7.
Suppose a lake is polluted with a certain chemical. If no more of this is added from some time onwards, how long will it take before the concentration drops down to 10 % of the value at that time ?
(To have a quantitative idea: for Lake Erie V ≈ 458 ((km )^3) and F ≈ 175 (km^3)/ year while for Lake Superior V ≈ 12221 ((km )^3)  and F ≈ 65,2 (km^3) year).

The tippler again.

The break down of alcohol in the liver is no match to the consumption when the latter is high, Consider

	(dx/dt)= v - (Subscript[g, out])(x)
	
where v =  (Subscript[g,  in]) and (Subscript[g, out ])is as defined as above:

	(Subscript[g, out])(x) = ((c x)/(k  +  x))

Exercise 8.
a) What does this model predict for  v<c  ?  Make a sketch of x(t).
b) And what for  v>c  ? Again make a sketch of x(t).
c) Do you agree that it is questionable in this last case whether the drinker will be able to keep v constant ?
d) Express both the steady state and the relaxation time in terms of the
parameters.