Archive

Posts Tagged ‘Calculus’

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 = []
sol.append(y0)
time = []
time.append(t)
#start of the program
while t<t_end:
	a = []
	a.append(sol[-1])
	for k in range(1,order):
		if k == 1:
			a.append(a[k-1]*(1-a[k-1]))
		else:
			sum = 0.0
			for j in range(0,k-1):
				sum = sum + a[j]*(-a[k-1-j])
			a.append((sum+a[k-1]*(1.0-a[0]))/k)
	sum = 0.0
	for i in range(0,order):
		sum = sum + a[i]*h**i
	sol.append(sum)
	t = t + h
	time.append(t)
#end of the program

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

#plotting the error versus time
sol = np.array(sol)
real_sol = np.array(real_sol)
time = np.array(time)
error = sol-real_sol
plt.plot(time[:],error[:])
plt.show()

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,

y(0)=0.

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.

A geometric proof of trigonometric derivatives

May 19, 2009 5 comments

As we already know, the first derivatives of \sin(x) and \cos(x) are \cos(x) and -\sin(x) respectively. The question is how do we prove it?

Since high school, the only proof that I knew is to use two of the trigonometric identities which are:

  • \sin \alpha - \sin \beta = 2 \sin (\frac{\alpha-\beta}{2}) \cos (\frac{\alpha+\beta}{2}), and
  • \cos \alpha - \cos \beta = -2 \sin (\frac{\alpha-\beta}{2}) \sin (\frac{\alpha+\beta}{2}).

Hence,

\frac{d}{dx}\sin x = \displaystyle{\lim_{\Delta x \rightarrow 0} \frac{\sin(x+\Delta x)-\sin(x)}{\Delta x} = \lim_{\Delta x \rightarrow 0} 2 \sin(\frac{\Delta x}{2})\cos(\frac{2x+\Delta x}{2})=\cos x}.

The derivative of \cos x can also be derived using the second identity.

However, there is another proof, which is (I think) more elegant. It is using a geometric interpretation of \sin x and \cos x. Let us draw a unit circle in \mathbb{R}^2.

geo_proof1

Thus we have y=\sin \theta and x=\cos \theta. We are looking for \frac{dy}{d \theta} and\frac{dx}{d \theta}, however we are not going to discuss the first derivative of \cos \theta as it can be done in the same way as \sin \theta.

\displaystyle{\frac{d \sin \theta}{d \theta}=\lim_{\Delta \theta \rightarrow 0} \frac{\sin(\theta + \Delta \theta)-\sin(\theta)}{\Delta \theta}=\lim_{\Delta \theta \rightarrow 0} \frac{y+\Delta y - y}{\Delta \theta}=\lim_{\Delta \theta \rightarrow 0} \frac{\Delta y}{\Delta \theta}}.

We are now going to estimate \Delta y, see the second figure.  geo_proof2The following relationship applies:

\sin \phi = \frac{\Delta y}{r} . or \Delta y = r \sin \phi.

Thus,

\displaystyle{\frac{d \sin \theta}{d \theta}=\lim_{\Delta \theta \rightarrow 0} \frac{\Delta y}{\Delta \theta}=\lim_{\Delta \theta \rightarrow 0} \sin \phi \frac{r}{\Delta \theta}},

where r is the chord PQ and \Delta \theta is the arc PQ. As \Delta \theta goes to zero, the ratio of the arc PQ and the chord PQ tends to 1 and the chord PQ tends to become a tangent line of the unit circle at (x,y), therefore \sin \phi = \sin (\theta + \pi/2). Assuming that the limit of multiplication is the same as the multiplication of limit, we have the following:

\frac{d \sin \theta}{d \theta}=\sin(\theta + \pi/2) = \cos(\theta).

As mentioned before, the first derivative of \cos \theta can be done the same way.

Also, I should thank Pak Yudi Soeharyadi as I got the idea of this geometric proof from him.

Fundamental Theorem of Calculus is fundamental !

May 11, 2009 1 comment

Most of my students do not seem to understand the two concepts of integrals. I have been there as well, so I did not blame them for not understanding the concept of integrals.  Even I sometimes forget. From that reason I try to post what the fundamental theorem of calculus means exactly.

It is all started with the definition first definition of integrals, which is formulated by the famous mathematician Bernhard Riemann. This is called the definite integral, or \int_a^b f(x) dx .

We can think that \int_a^b f(x) dx as the area below the curve f(x), see the picture below.

kalkulusThus,\int_a^b f(x) dx is the area S. As a first approximation, we shall divide the interval (a,b) into an equal number of segments, say \Delta x. Then we can approximate the area of S as follows:

\int_a^b f(x) dx \approx \sum_{i=0}^n f(x_i) \Delta x.

The right hand side of the equation above is called the Riemann sum. Our final step to find the value of definite integral is to take limit the Riemann sum as \Delta x \rightarrow 0,

\int_a^b f(x) dx = \lim_{\Delta x \rightarrow 0} \sum_{i=0}^n f(x_i) \Delta x.

Now, the second concept about integrals will be explained. This concept is called the indefinite integral or is also usually called the anti-derivative. Suppose F,f are functions of x such that F'(x)=f(x) (or we can say that the derivative of F(x) is f(x) or, the anti-derivative of f(x) is F(x)), then

\int f(x) dx = F(x) + C,

where C is a constant real number.

Hence, these two are two different concepts that are not related one to another. However, there is a great theorem that relates those two. The theorem perfectly and comfortably links the limit of Riemann sum and the anti-derivative. Actually there are two great theorems.

First Fundamental Theorem of Calculus Suppose f is a continuous function defined on closed interval [a,b] and define F as

F(x) = \int_a^x f(t) dt,

then F is continuous and differentiable, moreover, F'(x) = f(x).

Second Fundamental Theorem of Calculus Supposef is a function defined on closed interval [a,b] and F is the anti-derivative of f. If f is integrable on [a,b], then

\int_a^b f(x) dx = F(b) - F(a).

Therefore, if we want to compute, for examples, \int_0^1 x^2 dx, we don’t have to find the limit of Riemann sum of \sum_{i=0}^n x_i^2 \Delta x as \Delta x goes to zero.

Categories: Justincase, Teachings Tags: ,

Result of mid semester exam 2 of Calculus II (UTS 2)

May 7, 2009 2 comments

As I promised yesterday, you can now see the result of mid semester exam 2 of calculus II.

Please go to this page.

Next wednesday is gonna be our last tutorial. I hope you can do well in the final examination.

Categories: Teachings Tags: , ,

Result of mid calculus II mid semester exam (UTS 1)

March 18, 2009 10 comments

To all of my students in my tutorial class, subject Calculus IIa (T-19), Bu Mul and I finally finished the marking of your mid semester exams.

Klik here to go to my teaching page and go to Calculus IIa to find your exam results.

Categories: Teachings Tags: , ,