Posts Tagged ‘teaching’

Taylor method with automatic differentiation – 2

July 20, 2016 Leave a comment

In the previous post, we have talked about Taylor method to solve ordinary differential equation (ODE) numerically. All the derivatives in the Taylor method are computed using automatic differentiation. However, the drawback of this method is when we change our ODE, we have to write another code. In this post, we are going to present a way such that when we want to compute another ODE, we don’t need to change the majority of our code.

Continuing our analysis in the last post, we are going to need operations on polynomials. We are going to define operations such as:

  1. Addition between a constant and a polynomial
  2. Multiplication between  a constant and a polynomial
  3. Addition between polynomials
  4. Multiplication between polynomials
  5. Dividing a constant by a polynomial
  6. Dividing a polynomial by a polynomial
  7. Logarithm of a polynomial
  8. Exponentiation of a polynomial
  9. Power of a polynomial

Operation 1 and 2 are trivial. Suppose x(h) is a polynomial, then x(h)=\sum x_i h^i for some real constant x_i. Suppose \alpha is a real constant, then the addition between \alpha and x(h) is a polynomial y(h) where y_0 = x_0 + \alpha and y_i = x_i for all i \geq 1. The multiplication between \alpha and x(h) is a polynomial z(h) where z_i = \alpha x_i for all i.

Operation 3 and 4 are not difficult. Suppose x(h), y(h) are polynomials. Then the addition between x(h) and y(h) is a polynomial z(h) where z_i = x_i + y_i for all i. The multiplication between x(h) and y(h) is a polynomial w(h) where w_i = \sum_{j=0}^{i} x_j y_{i-j}.

Operation 5 and 6 are a little bit tricky. Suppose x(h) is a polynomial and \alpha is a constant. We want to compute y(h) where y(h) = \alpha / x(h). Notice that x(h)y(h)=\alpha. Therefore, we are going to take advantage of this equality to find the coefficient of y(h). Thus we shall have y_0 = \alpha / x_0, and y_i = - (\sum_{j=1}^{i} x_j y_{i-j})/x_0. For operation 6, let’s suppose z(h) = y(h) / x(h). As before, notice that x(h)z(h) = y(h), thus we can compute the coefficient of z(h) using this equality. Thus, z_0 = y_0 / x_0 and z_i = (y_i - \sum_{j=1}^{i} x_j z_{i-j} )/x_0.

Now, let’s consider the next operation, which is a logarithm of a polynomial. Suppose x(h) is a polynomial and we are going to find y(h)=\log(x(h)). Let’s differentiate both side with respect to h, then we have y'(h) = x'(h) / x(h). Thus by exploiting a division operation above, we have y_0 = \log x_0 and y_i = (ix_i - (\sum_{j=1}^{i-1} x_j ((i-j)y_{i-j})))/(ix_0).

Now, let’s consider the exponential operator, which is an exponentiation of a polynomial. Suppose x(h) is a polynomial and we are going to find y(h)=\exp(x(h)). Let’s take logarithm on both side and then differentiate it, we then have y'(h)/y(h) = x'(h), or y'(h)=y(h)x'(h). Therefore, by exploiting a multiplication operation, we shall have y_0  = \exp x_0 and y_i = (\sum_{j=0}^{i-1} y_j ((i-j)x_{i-j}))/i.

For the last operation of this post, we shall consider the power of a polynomial. Suppose x(h) is a polynomial, we are going to find y(h)=(x(h))^\alpha for any real number \alpha. First, we need to take logarithm on both sides and then differentiate the result with respect to h, y'(h)/y(h) = \alpha x'(h)/x(h) or y'(h) = \alpha x'(h)y(h)/x(h). By exploiting multiplication and division, we will get the coefficient of y(h).

To facilitate these operation, we build a Python module that will handle these operations. We will also write another code to integrate an ODE. This integration code will work in general without any dependencies on the system or on the order of the system.

def integrate(F,y_init,t_init = 0.0 ,tStop = 1.0 ,h = 0.05 ):
    dim = len(y_init)
    y_sol = []
    for i in range(0,dim):
    t_sol = []
    t = t_init

    while t < tStop:
        y = []
        for i in range(0,dim):

        for k in range(0,20):
            dy = F(y)
            for i in range(0,dim):
        n = len(y[0])
        for i in range(0,dim):
            temp = 0.0
            for j in range(0,n):
                temp = temp + y[i][j]*h**j
        t = t + h

    return np.array(t_sol) ,np.transpose(np.array(y_sol))

All the operations above and the last code will be written in a python file so that if we need to use this Taylor method, we can just import it. Let’s take an example to integrate one ODE. We will create another code that specifies the ODE and also imports all the operations above and then numerically integrate it.

from mytaylor import *
import matplotlib.pyplot as plt
import numpy as np
import math as m
#define your ODE here
def fun(y):
    z = subtractpoli(y[0],multiplypoli(y[0],y[0]))
    return [z]
#define your computation parameter here
y_init = [0.5]
t_init = 0
h = 0.1
tStop = 5.0
#integrating ODE
T, Y = integrate(fun,y_init,t_init,tStop,h)
#real solution computation
real_sol = []
for t in T:
	temp = m.exp(t)/(m.exp(t)+1.0)
#plotting the error versus time
real_sol = np.array(real_sol)
error = Y[:,0]-real_sol

In the above code, we try to numerically integrate y' = y(1-y) with initial condition y(0)=0.5. This is the same system we try to integrate in the previous post. I hope you can see the difference between the last code and the code we have shown in the previous post. We really don’t need to consider about Taylor method at all, we just need to write our system in the above code if we want to compute the numerical solution of our code. Notice that we have written all the operations and the integrating code in a file named which is imported in the first line of the above code. In the fifth line, we can write any ODE we want to integrate, however, the way we write our ODE will be slightly different. Due to polynomial operations that we have defined in this post, x+y cannot be written as just x+y but addpoly(x,y). Once we define our ODE, we can just integrate it like in line 15.

Another thing we would like to do is to extend the operation of polynomial. We have not included trigonometry operations on polynomials and et cetera. We would like also to test this Taylor method to a stiff system. Please find the following two files ( and for interested readers if you would like to test yourselves.


Taylor method with automatic differentiation

June 14, 2016 1 comment

I teach numerical analysis course and in the one of the chapters, I teach my students about the numerical solution of ordinary differential equations (ODEs). One of numerical solution methods I teach is the so-called Taylor method. This method works by applying the Taylor series to the solution of the differential equation. Suppose we have one-dimensional differential equation as follow:

\frac{dy}{dt} = f(y) and y(t) = y_0.

Let’s assume that we only consider an autonomous ODE first.  To find a solution at the next step, y(t+h), we need to apply Taylor expansion, as follow:

y(t+h) = y(t) + y'(t) h + y''(t)\frac{h^2}{2} +y'''(t)\frac{h^3}{6}+\dots.

Therefore, depending on how accurate we want our solution is, all we need is to compute the derivatives of y(t). Unfortunately, computing all the derivatives is not an easy task for a computer. We have to derive all the derivatives by ourselves and then input them to the computer.

In this note, I will introduce an automatic differentiation as a tool to compute the derivatives of y(t+h) so that we can compute y(t+h) as accurately as possible. Let us consider the series expansion of y(t+h) as above, but I will use constant a_i to replace the i-th derivative of y(t),

y(t+h) = a_0 + a_1 h + a_2 h^2 + \dots.

Let us derive the above series with respect to h.

y'(t+h) = a_1 + 2 a_2 h + 3a_3 h^2 + \dots.

Notice that the derivative at the left hand side of the above equation is actually equal to f(y(t+h)) since we know that y(t+h) satisfies the ordinary differential equation. Then we have,

f(y(t+h)) = a_1 + 2 a_2 h + 3a_3 h^2 + \dots.

Since y(t+h) is a series of h then f(y(t+h)) is also a series of h. Then both sides of the above equation are series of h, thus we can equate each of the coefficient of h^j to get a_1a_2 and so on. The most difficult task in this problem is to find the series of f(y(t+h)) since the function f(y) can be anything.

Let’s take an example. Suppose we want to integrate the following differential equation,

\frac{dy}{dt} = y(1-y) and y(0)=0.5.

Let us choose h=0.1 and compute y(0.1). As before, we assume that y(0.1) is a series of h as follows,

y(0.1) = a_0 + a_1 h + a_2 h^2 + \dots.

Differentiating both sides of the above equation, we get

y'(0.1) = a_1 + 2 a_2 h + 3 a_3 h^2 + \dots.

The right hand side of the above is equal to

y'(0.1) = y(0.1)(1-y(0.1)) = (a_0 + a_1 h + \dots)(1  - (a_0 + a_1 h + \dots)) = b_0 + b_1 h + b_2 h^2 + \dots.

Therefore, we get,

a_1 = b_0a_2 = b_1/2,a_3 = b_2/3,  and so on.

The real question is how to find b_0, b_1 and so on. If the function f(y) only involves operations such as addition, subtraction and multiplication, it would be easy to find the b‘s. But if the function f(y) has operations such as division, exponential, logarithm and trigonometric functions then we have to think a bit harder on how to find b‘s. We go back to the previous example. Since the operations in this example are only multiplication and subtraction, we can easily compute b_0, b_1 and so on. Using algebraic modification, we get the following,

b_0 = a_0 \times (1-a_0), thus a_1 =a_0 \times (1-a_0)
b_1 = (a_0 \times -a_1 ) + (a_1 \times (1-a_0)), thus a_2 =((a_0 \times -a_1 ) + (a_1 \times (1-a_0)))/2
b_2 = (a_0 \times -a_2 ) + (a_1 \times (-a_1))+(a_2 \times (1-a_0)), thus a_3 =((a_0 \times -a_2 ) + (a_1 \times (-a_1))+(a_2 \times (1-a_0)))/3

and so on.

The more a‘s we compute, the more accurate our solution is. Using Taylor method this way, we can easily apply this method in a program. It should be clear why this method is called automatic differentiation because we don’t need to manually derive all the derivatives.

Next, I tried to numerically integrate the previous example. Below is the graphic of the error of Taylor method with automatic differentiation versus time. As we can see from the figure, that the error of our numerical integrator is so small and it is almost exact. The number of a I used in this graph is 10. error_plot

I use python to numerically integrate this example. You can find the code to produce the above figure in the end of this note. One of the drawbacks of this method is when you change your ODE, you have to change your code as well. In other numerical integrators such as Runge Kutta, we don’t need to change the whole code, we just need to change the ODE and we should get the solution. Together with my student, I am trying to tackle this problem and I hope it will be posted in here.

import math as m
import numpy as np
import pylab as plt

#preliminary constants
order = 10
h = 0.1
t = 0.0
t_end = 5
y0 = 0.5
#defining solution
sol = []
time = []
#start of the program
while t<t_end:
	a = []
	for k in range(1,order):
		if k == 1:
			sum = 0.0
			for j in range(0,k-1):
				sum = sum + a[j]*(-a[k-1-j])
	sum = 0.0
	for i in range(0,order):
		sum = sum + a[i]*h**i
	t = t + h
#end of the program

#real solution computation
real_sol = []
for t in time:
	temp = m.exp(t)/(m.exp(t)+1.0)

#plotting the error versus time
sol = np.array(sol)
real_sol = np.array(real_sol)
time = np.array(time)
error = sol-real_sol

How to fold a paper into thirds, a real analysis perspective

February 19, 2016 Leave a comment

Folding a paper into third is useful when you want to mail someone using an envelope. I search around the internet and there are so many ways to do this. So far, I didn’t see the method I am about to explain. I found this method by accident while I am teaching Real Analysis this semester. I am going to prove that this method actually works. Other methods that I found always use Geometry in the proof but I am going to use Real Analysis instead.

The procedure

Please have a look at the following figure. firststep

The above figure is the first step of our method. Take a sheet of paper and let’s call the left side, side A and the right side, side B. First, you need to fold the paper side A as in the second image of the figure above. It doesn’t have to divide the paper equally, you can just fold it anywhere. Once you do that, you undo the fold and you will have a fold mark, which we call it x_0. Next, you need to fold again, but now you need to fold the right side (side B) along the fold mark x_0 as in the second row of the above figure. Now you will have two fold marks, the latter is called y_0. Now take a look at the following figure. secondstep

Now that we have two fold marks (x_0 and y_0), we need to fold side A along the fold mark y_0, as in the first row of the above figure. Undo the fold, you will have the third fold mark (we call it x_1). Next, we fold the side B until it touch the fold mark x_1 and we will have four fold marks, x_0, x_1 and y_0, y_1. Continue the step over and over again and at one stage, at the n-step, x_n will be close to one third of the paper and y_n will be very close to two third of the paper.


Let’s take an example. Suppose in the beginning, let’s fold side A such that we divide the paper equally. Assume the length of a paper is 1. This means x_0=0.5. If we simulate this case, we will get the following table.

n 1 2 3 4 5
x_n 0.5 0.375 0.3475 0.3359375 0.333984375
y_n 0.75 0.6875 0.671875 0.66796875 0.666992188
6 7 8
0.333496094 0.333374023 0.333343506
0.666748047 0.666687012 0.666671753

As we see, that it took approximately 4-8 steps until x_n is closed to one third of a paper or y_n is closed enough to two third of the paper. However, what happens if we choose different initial condition. We simulate this using a number of initial condition. Suppose we use an A4 paper which has length of 29.7 cm. We choose a number of different initial conditions. We can see the results in the following figures. simulation01The x-axis from the above figure is the initial condition and the y-axis is the number of iteration needed until it is closed enough to one third of paper. Since we use an A4 paper, the x-axis ranges from 0 to 29.7 cm.  The tolerance we use in the above figure is 1 mm. In conclusion, we only need to do only maximum four steps until the last fold mark is only 0.1 cm away from one third of the A4 paper. simulation05

We relax the tolerance in the above figure. We use 0.5 cm as the tolerance and the number of iterations is slightly improved. It can be seen that if our first trial is close enough to one third of A4 paper, we only need to do one step.

Mathematical proof

Finally, it comes to prove that this methods should work. mathematically not numerically. First of all, we need to mathematically model the above procedures. Let’s assume that the length of the paper we use is 1. Let’s use the same notation as used in the above procedures. Given x_0 we will get y_0 = (1+x_0)/2. Recursively, we could compute x_n and y_n as follow:

x_n = y_{n-1}/2 and y_n=(1+x_n)/2, as n=1,2, … .

Even though there are two sequences (x_n) and (y_n), we only need to consider (x_n). We could simplify the sequence (x_n) as follow.

x_0 is given and x_n=y_{n-1}/2=((1+x_{n-1})/2)/2=(1+x_{n-1})/4 for n=1,2,\dots.

Lemma 1
When 0<x_0<1/3 (as 1/3<x_0<1), the sequence (x_n) is bounded above (bounded below).

Lemma 2
When 0<x_0<1/3 (as 1/3<x_0<1), the sequence (x_n) is monotonically increasing (decreasing).

Theorem 3
Given 0<x_0<1, the sequence (x_n) converges to 1/3.

Corollary 4
Given 0<x_0<1, the sequence (y_n) converges to 2/3.

The proofs of the above results are given to my students as exercises.

This is the end of this note and the following is the code to produce the above figures.

import math
import numpy as np
import pylab as plt

def f(x):
 return (29.7+x)/2.0

def g(y):
 return y/2.0

def number_of_iter(x,tol,max_iter):
 iter = 0
 er = abs(x-29.7/3)
 if er <= tol: return 1 
 while er > tol and iter < max_iter: 
  iter = iter + 1 
  y = f(x)
  x = g(y)
  er = min(abs(29.7/3-x),abs(2*29.7/3-y))

 return iter

#main code
x = np.arange(0,29.7,0.1)
y = []
tol = 0.5
max_iter = 100

for t in x:
 result = number_of_iter(t,tol,max_iter) 
y = np.array(y)
ax = plt.gca()

How to compute distance between a point and a line

January 27, 2016 Leave a comment

This semester (Even Semester, 2015-2016) I am teaching a new course, that is called Analytic Geometry or Geometry Analytic, both are the same, I think. Don’t get me wrong, the course is not new, it is just I would be teaching this course the first time. When they told me that I was going to be the lecturer of this course, they did not give me a standard, they did not give me a textbook they normally use, or a list of topics that I must cover, perhaps because this course is not a compulsory course. So, I have a freedom, I can choose topics that I want to teach in Analytic Geometry. The first thing I did, was to browse the internet any textbook about Analytic Geometry or any lecture note. In the end, I pick a textbook (Indonesian book) and I choose some topics of my interest to be the material of this course.

Okay, enough for the background, I am sure you don’t want to hear anymore on that. Let’s get back to main point of this note. Long story short, a straight line became one of the topics in my course, and I was interested on how the formula to compute distance between a point and a line is derived.In this case, I only consider two dimensional case and by distance, I mean the shortest distance between such a point and a certain line, which also means the perpendicular distance. If you don’t remember what is the formula, let me recall you.

Consider a straight line l, given by the following equation:

l : ax + by + c=0.

Suppose we have a point P(x_p,y_p), then the distance between the point P and the straight line l is

d = \frac{|ax_p+by_p+c|}{\sqrt{a^2+b^2}}.fig1-crop

There are other proofs on how to derive the above formula, but I really really like the proof I am going to show you below.  I will divide this note into two section. The first section will talk about a straight line and how to get an equation of a straight line. The latter section will derive the formula.

1. A straight line equation

In high school, I think you must know what conditions we need to have in order to obtain a straight line equation. If we have two points in XY-coordinate we can compute the line equation using the following:


If we have a gradient and a point, we can also compute the line equation using the following:


Of course, in a problem, we won’t get this information so easily. We have work a bit more to get either two points or one point and a gradient and then we are able to obtain the equation of the straight line.

Before I teach this course, I can only conclude that if you want to know the equation of a straight line you need to know either:
(i) two points, or
(ii) a point and gradient of the line.
But, actually there is a third condition, and if we know this condition, we can also compute the equation of a straight line. In the third condition, a line is assumed to be the tangent line of a certain circle. Thus, the characteristics we need to know to form a line are the radius of the circle and the angle between the radius and the positif x-axis. See the figure below.fig2-crop

In the figure, we can see that a line can be determined uniquely if the radius of the circle n and the angle \alpha are known. In the first section of this note, we shall derive how to define a line equation given these two conditions.
Suppose we have a straight line, n are \alpha are given. See the figure below. Consider the point P(x,y). We are going to find the relationship of x and yfig3-crop

The radius n, which is |ON|,  can be computed by adding |OM| and |MN|. We shall consider |OM| first. Consider triangle OMQ which is a right triangle at M. We have the relationship:

|OM|=|OQ|\cos \alpha=x \cos \alpha.

Consider another triangle, PM'Q, which is also a right triangle at M'. We have the following relationship

|PM'|=|PQ|\sin \alpha=y \sin \alpha as it can be computed that the angle \angle PQM is also \alpha.

Therefore, we have n=x \cos \alpha + y \sin \alpha, which is the equation of straight line given that n and \alpha is known.

Before we discuss how to derive the formula to compute the distance between a point and a line, I would like to discuss how to find n and \alpha if we know the general equation of a straight line ax+by+c=0. Consider a line equationax+by+c=0, then we move $c$ to the right hand side and multiple both sides with a non zero constant kkax+kby=-kc. We shall choose k such that ka=\cos \alpha and kb=\sin \alpha. Therefore, we have

k = \pm \frac{1}{\sqrt{a^2+b^2}} and

our equation of line becomes:


Here, we have that the right hand side of the above equation is n=\pm(-\frac{c}{\sqrt{a^2+b^2}}), and we choose the positive value to get n. The angle \alpha can also be computed once we determine n.

2. Distance between a point and a line

To compute a perpendicular distance between a point and a line, we shall use the above result. Consider a straight line l:x \cos \alpha+y \sin \alpha = n_1 and a point P in the following figure.fig4-crop


We want to compute d. In order to do that, we need to make another line that is parallel to the line l and passing through point P. This line, because it is parallel to l and, has the following equation:

x \cos \alpha+y \sin \alpha = n_2

Thus the distance between the point P and the line l can be easily computed by considering the absolute difference between n_1 and n_2 as follow:

d = |n_2-n_1|=|x_p \cos \alpha + y_p \sin \alpha - n_1|

As we know from the first section that we can substitute the latter expression into:

d = |x_p \frac{a}{\sqrt{a^2+b^2}}+ y_p\frac{b}{\sqrt{a^2+b^2}}-(-\frac{c}{\sqrt{a^2+b^2}} )|



which is the same as the formula we mentioned in the beginning of our note.

Categories: Teachings Tags: , ,

Notes on diagonalizable linear operators

May 11, 2011 3 comments

Okay, I have a few things in my mind about this topic and I don’t want to lose them tomorrow or the day after. Therefore, I created this note. Hope it is useful to other people as well. Background of this note are as follows:

  1. Suppose A is an n \times n matrix and it is diagonalizable. Therefore, by definition, there is an invertible matrix P and a diagonal matrix D. However, I always forget the relationship between A and D, whether it is D=PAP^{-1} or D=P^{-1}AP. Generally, those two are no different, but we can choose P such that the column vectors of matrix P are the eigenvectors of matrix A, hence we have to know precisely which one it is, otherwise it gets wrong. Until now, when I forget about this, I always derive it and it takes some of my precious time (haha) but now after I have this, I can just open the Internet and look for this note.
  2. Point 1 above talks about diagonalizability of a matrix. But generally, we have a linear operator acting on a general vector field instead. A linear operator is called diagonalizable if there is a basis B such that the matrix representation of this linear operator is a diagonal matrix. It seems for me, diagonalizability of a linear operator and that of a matrix are talking about two different things. One is talking about how to find a basis, while the other is talking about how to find an invertible matrix. However, these two concepts turn out to be the same.

We begin with the following theorem.

Theorem 1 (Coordinate-change matrix)
Let E be a standard basis and B be another basis for an n-dimensional linear space V. Then there is an invertible matrix P such that for every v \in V

P (v)_B=(v)_E.                                                                                                          (1)

Matrix P is called the transition or coordinate-change matrix from the basis B to the basis E. Note that (v)_B and (v)_E are the coordinates of vector v with respect to basis B and E respectively. In fact, the column vectors of matrix P are the coordinates of the basis vectors of B with respect to the standard bases E (i.e. if B=\{ v_1,v_2,\dots,v_n \} then P=[(v_1)_E (v_2)_E \dots (v_n)_E] ).

We won’t prove the theorem as it can be found in any linear algebra textbook. However, using this theorem we can relate the diagonalizability of a linear operator and that of a matrix.Suppose T is a diagonalizable linear operator on an n-dimensional linear space V. Take x \in V then let y be y=T(x). Suppose E is a standard basis in V and B is a basis in V that consists of independent eigenvectors of T.

Then we have the following.

(i) (y)_E = [T(x)]_E = [T]_E (x)_E = A (x)_E, where A is a representation matrix of T with respect to standard basis E.
(ii) (y)_B = [T(x)]_B = [T]_B (x)_B = D (x)_B, where D is a representation matrix of T with respect to basis B. We also know that D is a diagonal matrix.

We want to show that using the theorem 1, D=PAP^{-1} or D=P^{-1}AP. From equation (1) we get:

P (y)_B = (y)_E = A (x)_E = A P (x)_B.

Hence, we have (y)_B = P^{-1}AP (x)_B. While from point (ii) we have (y)_B = D (x)_B. As a conclusion, we have,

D = P^{-1} A P.

We have derived that in fact, D = P^{-1} A P is the right one. We also derived that the invertible matrix P has such a close relationship with the basis B as the column vectors of P are the basis vectors of B, which are the eigenvectors of the linear operator T. This conclude my note.

Academic Recharging

May 6, 2011 Leave a comment

I first hear about this phrase in the Indonesian government of education website. Academic recharging is one of fund application schemes provided to academician who feels bored and sometimes tired doing teaching in their respective university. This scheme is provided to ‘recharge’ them in some way such that they are able to go abroad conducting research with their colleagues overseas. Indonesian government hopes by doing this scheme, lecturer and academician could update their knowledge to the current topics so that Indonesian academic society would not be left behind.

But I have found my way of academic recharging. And I give many thanks to ICTP and MIT. Courses that are provided by ICTP and MIT have been recorded and have been available online. I don’t know since when this happens but for me it is great. I can update my knowledge with some new information that I have never learned before.

After downloaded a video, I can watch the course that I am interested. When I don’t understand a thing, I can just pause and google them. I can pause anytime I want when some students knocking my door. I strongly recommended this way of learning to my colleagues and my students as this is much as fun as watching movie.

Here are some references that you might like.

  1. http://

Mixture problem in Calculus

April 29, 2010 7 comments

I keep forgetting the way to solve this problem. So, when I in future have no clue I know where to look.

The problem
A 200-galon tank is half full of distilled water. At time t=0, a solution containing 0.5 pound/galon of salt enters the tank at the rate of 5 gal/minutes, and the well-stirred mixture is withdrawn at the rate of 3 gal/minutes.
a. Give a mathematical model representation of the above problem.
b. At what time, will the tank be full?
c. At the time the tank is full, how many pounds of salt will it contain?

The mathematical model
The model describing the above process is based on the following formula:

rate of change of the amount of salt  = (rate of salt arrived) – (rate of salt departed)

Suppose V(t) is the volume of the liquid in the tank at time t. The volume of liquid inside the tank will increase as the inflow (5 gal/minutes) is bigger than the outflow (3 galon/minutes).  Thus, we have

V(t)= 100 galon + (5 gal/min – 3 gal/min) \times t min,
V(t)=(100 + 2t) gal,

as we know that initially the tank is half full (100 gal).

Suppose y(t) is the amount of salt inside the tank at time t. The rate in of salt is as follows,

rate in (pound/min) = (0.5 pound/gal) \times (5 gal/min) = 2.5 pound/min,

while the rate out is
rate out (pound/min) = \frac{y(t)}{V(t)} \times (pound/gal) outflow (gal/min)
= \frac{y}{100+2t} \times 3 (pound/min)
= \frac{3y}{100+2t} (pound/min).

Therefore, the mathematica model is given by,

\frac{dy}{dt}=2.5 - \frac{3y}{100+2t}.          (1)

The solution
First we determine time t at which the tank will be full. We use the following equation,

V(t)=(100 + 2t),

We solve for t when the volume V(t) is 200 galon,  then we obtain t=50.

To solve the differential equation (1) we need the initial condition
as follow,


We are going to use Maple to solve the problem

> model := diff(y(t),t) = 2.5 - (3*y(t))/(100+2*t);
> init := y(0)=0;
> dsolve({model,init});

We obtain the solution that is given by

y(t) = 50+t-\frac{2500 \sqrt{50}}{(50+t)^{3/2}}

When t=50, the amount of salt inside the tank will be

y(50) = 82.32233047.