My *nix world
home > Maths > Calculus > Solving ordinary differential equations (ODE)

Solving ordinary differential equations (ODE)

Introduction

In Algebra we have what is called algebraic equation that is an mathematical expression described in terms of an unknown variable x.

For instance 2x+7=0 is an algebraic equation of first order(because the highest order of the x term is 1) while x^5-3x+1=0 is an algebraic equation of fifth order (because the highest order of the x term is 5). The solution of an algebraic equation (if it really exists) belongs to the real set or to the complex set of numbers. Note that due to the fact that \mathbb{R}\subset\mathbb{C} we could simply say that the solution (if exists) must be a complex number (which in some particular cases could be a real number, a rational number, an whole number or an natural number).

So an algebraic equation is expressed in terms of x-variable and to find its solution(s) means to find a number for that x-variable such that it verifies the equation (when it is used the left expression is equal with the right expression).

Very similar to an algebraic equation a differential equation(DE) is also an equation that is expressed in terms of y-function (instead of x-variable) and thus the solution of a differential equation is not a number (like in the case of the algebraic equation) but a function. Moreover, a differential equation is not only expressed in terms of the y-function but also (and this is important) in terms of its derivative, i.e. \frac{dy}{dx}. That's why it is called differential in the first place, it must include the derivative. The higher the order of derivative the higher is the order of differential equation. For instance, y''+2y'-y=0 is a DE of the 2nd order (the higher order of derivative of y is 2). Note that such equation can be expressed in the same time in terms of x. For instance y''+2y'-y=3x+7 is also a DE.

Note that while an algebraic equation might have only a solution that verifies the equation a DE might have an infinity of solutions that verify the DE. Why is that? Let's imagine the function that is associated to a algebraic equation. The solution of that equation are those x-values where the function cross the x-axis. Now let's image a differential equation. A solution of that DE is a function (whatever that means). If we shift the function (the graph of the function) just a bit to the left or to the right, or even a bit higher or lower we still have the same function just that it was a bit shifted away (from the x,y-axes). Each shifted version of that function is a (particular) solution of the DE. So by drawing the graph of the function and by shifting around we get a multitude of solutions (an infinity) although the function graph will still look the same.

An ordinary differential equation (ODE) is an differential equation expressed only in terms of only one variable (like f(x), we express f in terms of one variable x). If it is expressed in terms of many variables (like f(x1,x2,...)) then it is called partial differential equation (PDE) and that's something else which I am not going to develop here.

When talking about equations (especially about ODE) sometimes we hear "homogeneous" and "non-homogeneous" equations. "Homos genos" comes from Greek and means "of the same kind". Thus being non-homogeneous means they belongs to a distinct kind, like functions (y) and variables (x).

Now that you know what a differential equation is it would be also helpful in showing different methods that occurs in solving ordinary differential equations (ODE).

ODE of 1st order

So an ordinary differential equation where the highest derivative has the first order (i.e. y') is called ODE of 1st order.

The general form of a homogeneous ODE of 1st order is:

y'+a \cdot y=0. This ODE is expressed only in terms of function y which is something of the same kind, thus homogeneous.

We could rewrite that equation like this: y'=-ay. The question is: what should be that function y such that when is derived we get a new function (y') that looks exactly like the first one (i.e. y) except that it is multiplied by a constant (eg. -a) ? Of course, it is e^{-a\cdot x}. A good guess would be that our solution looks like this: y=e^{r \cdot x}.

To verify if our guess is valid would should find out its derivative and then we should use these in the original ODE Let's verify if the left term (LT) is equal with the right term (RT), i.e. LT=RT.

RT=0

LT=y'+a\cdot y=\frac{dy}{dx}(e^{-a \cdot x})+a \cdot (e^{-a \cdot x})=-a\cdot e^{-a\cdot x}++a \cdot e^{-a \cdot x}=0=RT

Because LT=RT we can conclude that our guess is valid, i.e. y=e^{-a \cdot x} is a solution of the ODE. But this is just one solution because if we guess that our solution would be 100\cdot e^{-a\cdot x} it will be also valid. So the general solution of such homogeneous ODE is:

y=C\cdot e^{r\cdot x} where C is an arbitrary constant and -r is the coefficient of the y-term of the ODE.

The general form of a non-homogeneous ODE of 1st order is:

y'+a \cdot y=f(x) where f is any function of x (eg. f(x)=3x+1). Because this equation mixes both functions (y) and variables (x) such that we can distinguish between them it is said to be a non-homogeneous ODE.

If we assume that the solution for such an equation is a function y=C\cdot e^{r\cdot x} it would be wrong because it cannot be verified. To solve such an equation means to find the solution y_h of the homogeneous part of the ODE and also to find the solution y_p of the f(x) part of the ODE. Thus the general solution of the ODE would be y=y_h+y_p

Example:

y'-2y=1+x

Note that the same equation could be written as y'=2y+1+x or \frac{dy}{dx}=sy+1+x. It is just the same equation, maybe written using a different notation. That's all.

Step 1: find the solution y_h of the homogeneous part of the ODE: y'-2y=0.

We suppose that y_h=C\cdot e^{r\cdot x} is the general solution of the homogeneous equation.Its derivative is y_h'=C\cdot r\cdot e^{r\cdot x}. If we use these in the homogeneous part of ODE (i.e. y_h'-2y_h=0) we get:

C\cdot r\cdot e^{r\cdot x}-2 \cdot (C\cdot e^{r\cdot x})=0 <=> C\cdot e^{r\cdot x}(r-2)=0. The product is null if any of these terms are null, i.e. C\cdot e^{r\cdot x}=0 or r-2=0. This second equation is also called "the characteristic equation of ODE". If we assume that C\neq 0 then C\cdot e^{r\cdot x}\neq 0\forall x,r,C. So r-2 must be noll, i.e. r=2. Thus y_h=C\cdot e^{2x} is the general solution for the homogeneous part of the ODE.

Step 2: find the solution y_p of f(x) part of the ODE

As we see f(x) is the first order algebraic equation and thus it can be written like this: y_p=ax+b where a and b are two arbitrary constants that we are going to determine right now.

y_p'=a. If we use these in the f(x) part of the ODE (i.e. y_p'-2y_p=1+x) we get:

a-2\cdot (ax+b)=1+x <=> (a-2b)-2a\cdot x=1+x. We equalize the non x-terms from both sides and likewise the x-terms from both sides and solve that system of equations:

\begin{cases}(a-2b)=1\\-2a=1\end{cases} => a=-0.5 and b=-0.75. Also y_p=-0.5x-0.75

Step 3: the general solution of the ODE is y=y_h+y_p

y=C\cdot e^{2x}+(-0.5x-0.75)=C\cdot e^{2x}-0.5x-0.75 where C is an arbitrary constant. For a particular value of C you'll get a particular solution of y.

ODE of 2nd order

So an ordinary differential equation where the highest derivative has the second order (i.e. y'') is called ODE of 2st order.

Such an ODE if homogeneous looks like y''+ay'+by=0 or if it's non-homogeneous then it looks like y''+ay'+by=f(x). The way we solve this equation is similar to that we solve the ODE of 1st order. Its characteristic equation will be r^2+a\cdot r+b=0 which has the following solutions: r_{1,2}=\frac{-a\pm\sqrt{\Delta}}{2} where \Delta=a^2-4\cdot 1 \cdot b.

Depending on the value of \Delta the solution for our ODE can vary.

2 real solutions

When \Delta>0 the characteristic equation has two distinct real solutions, namely r_1 and r_2. Thus the homogeneous ODE has two distinct solutions:

y_{h1}=C_{1}\cdot e^{r_1 \cdot x} and respective y_{h2}=C_{2}\cdot e^{r_2 \cdot x}. Also the solution of our ODE is the sum of these two solutions, namely:

y_h=y_{h1}+y_{h2}=C_{1}\cdot e^{r_1 \cdot x}+C_{2}\cdot e^{r_2 \cdot x}.

1 double solution

When \Delta=0 the characteristic equation has one double real solution, namely r_1=r_2=r=-\frac{a}{2}. Thus the homogeneous ODE has the following solutions:

y_{h1}=C_{1}\cdot e^{r_1 \cdot x} and respective y_{h2}=C_{2}\cdot e^{r_1 \cdot x}. Also the solution of our ODE is the sum of these two solutions, namely:

y_h=y_{h1}+y_{h2}=C_{1}\cdot e^{r \cdot x}+C_{2}\cdot e^{r \cdot x}=e^{r \cdot x}\cdot(C_{1}+C_{2})=e^{r \cdot x}\cdot(C) where C=C1+C2 is just an arbitrary constant. But this is 2nd order ODE and thus it MUST have two solutions not just one. One solution is y_{h1}=C\cdot e^{r\cdot x} but what is the other solution?

Well, we assume that the second solution y_{h2} is similar to the first one y_{h1} multiplied by an unknown function v(x), i.e. y_{h2}=y_{h1}\cdot v(x)=C\cdot e^{r\cdot x}\cdot v(x).

y_{h2}'=y_{h1}'\cdot v +y_{h1}\cdot v'=C\cdot r\cdot e^{r\cdot x}\cdot v+C\cdot e^{r\cdot x}\cdot v'

y_{h2}''=C\cdot r^2\cdot e^{r\cdot x}\cdot v+C\cdot r\cdot e^{r\cdot x}\cdot v'+C\cdot r\cdot e^{r\cdot x}\cdot v'+C\cdot e^{r\cdot x}\cdot v''

If we use these two in the original ODE (i.e. y_{h2}''+a\cdot y_{h2}'+b\cdot y_{h2}=0) we get:

C\cdot r^2\cdot e^{r\cdot x}\cdot v+C\cdot r\cdot e^{r\cdot x}\cdot v'+C\cdot r\cdot e^{r\cdot x}\cdot v'+C\cdot e^{r\cdot x}\cdot v''+a\cdot(C\cdot r\cdot e^{r\cdot x}\cdot v+C\cdot e^{r\cdot x}\cdot v')+b\cdot(C\cdot e^{r\cdot x}\cdot v)=0 <=>C\cdot e^{r\cdot x}\cdot v''

This product is null only if v''=0. If we take integral of both sides we get that v'=C_1 where C1 is an arbitrary constant. We take again integral of both sides and we get v=C_1\cdot x+C_2 where C_2 is another arbitrary constant.

So our second solution is y_{h2}=C\cdot e^{r\cdot x}\cdot v(x)=C\cdot e^{r\cdot x}\cdot (C_{1}\cdot x+C_{2}). We assume that there exist two arbitrary constant A, B such that A=C*C1 and B=C*C2, thus y_{h2}=(A\cdot x+B)\cdot e^{r\cdot x} which is also the general case for y_{h1}.

So the general solution is y=y_{h2}=(A\cdot x+B)\cdot e^{r\cdot x}

2 complex solutions

When \Delta<0 the characteristic equation has two distinct complex solutions, namely r_1 and r_2. Thus the homogeneous ODE has two distinct solutions:

y_{h1}=C_{1}\cdot e^{r_1 \cdot x} and respective y_{h2}=C_{2}\cdot e^{r_2 \cdot x}. Also the solution of our ODE is the sum of these two solutions, namely:

y_h=y_{h1}+y_{h2}=C_{1}\cdot e^{r_1 \cdot x}+C_{2}\cdot e^{r_2 \cdot x}. But this solution is expressed in terms of complex numbers and this we don't like (it is not necessarily wrong but it would be cumbersome if we want to use this form to plot a graph). To reduce this form to something expressed only in terms of real number we could use the Euler's identity:

e^{i\cdot ax}=cos(ax)+i\cdot sin(ax)

So if the characteristic equation has the complex solutions: r_{1,2}=m\pm i\cdot n the ODE's solutions would be: y_h=C_{1}\cdot e^{(m+i\cdot n) \cdot x}+C_{2}\cdot e^{(m-i\cdot n) \cdot x}=C_{1}\cdot e^m\cdot e^{i\cdot nx}+C_{2}\cdot e^m\cdot e^{-i\cdot nx}=e^m\cdot cos(nx)\cdot(C_{1}+C_{2})+i\cdot e^m\cdot sin(x)\cdot(C_{1}-C_{2})

We note C1=p+iq and C2 its conjugate i.e. C2=p-iq. Thus y_h=e^m\cdot cos(nx)\cdot(2p)+i\cdot e^m\cdot sin(x)\cdot(2q\cdot i)=e^m\cdot cos(nx)\cdot(2p)+e^m\cdot sin(nx)\cdot(-2q). If we note A=2p and B=-2q then our solutions becomes:

y_h=e^m\cdot (A\cdot cos(nx)+B\cdot sin(nx)) where A,B are two arbitrary constants.

Numeric methods for ODEs

Sometimes we deal with such ODE that are difficult to solve if not impossible. In such situations we can use numerical methods to determine an approximative solution for ODE. There are many numerical methods to use but one I would prefer (due to its simplicity) is the Backward Euler method.

Backward Euler method

This method assume that:

(i) the ODE is 1st ODE

and

(ii) we know ahead a point that lies on the graph of our solution.

What we do is to use this starting point and then find out the value of the slope/tangent to the graph at that point. Then we move an increment bit (a step) in one direction (either the positive direction of x-axis or the negative direction of x-axis, depending which part of the graph we want to depict) by using the slope we've already computed and thus we get a new point. We repeat this procedure as long as we want. When we finished the procedure we've got a visual graph of how the ODE's solution looks like.

Note: the smaller the increment (the step) the better the function approximation.

Example:

y'=1+3x-2y given that the (1,2) is on the graph. We assume a step of h=0.25

Step 1: we mark our starting point (1,2) and then find out the slope k=y'. When x=1 and y=2 the slope k=y'=1+3*1-2*2=1+3-4=0. So the graph slope at (1,2) is parallel to the x-axis. We draw a small segment with the length h and slope k.

Step 2: we compute the new point coordinates:

x=1+h=1+0.25=1.25 and y=kh+2=0*h+2=2

Step 3: the new point became (1.25, 2) and the slope at this point is k=y'. When x=1.25 and y=2 the slope k=y'=1+3*1.25-2*2=0.75. We draw a small segment with the length h and slope k.

Step 4: we compute the new point coordinates:

x=1.25+h=1.25+0.25=1.5 and y=kh+2=0.75*h+2=2.1875

Step 5: the new point became (1.5, 2.1875) and the slope at this point is k=y'. When x=1.5 and y=2.1875 the slope k=y'=1+3*1.5-2*2.1875=1.125. We draw a small segment with the length h and slope k.

Step 6: we compute the new point coordinates:

x=1.5+h=1.5+0.25=1.75 and y=kh+2=1.125*h+2.1875=2.46875

Step 7: the new point became (1.75, 2.46875) and the slope at this point is k=y'. When x=1.75 and y=2.46875 the slope k=y'=1+3*1.75-2*2.46875=1.3125. We draw a small segment with the length h and slope k.

Step 8: we compute the new point coordinates:

x=1.75+h=1.75+0.25=2 and y=kh+2=1.3125*h+2.46875=2.796875

The graph drawn by these points looks like this:

backward-Euler-method

Using these points we could use a regression tool (like this) to determine the best fitting line. The greater is the determination coefficient (r) the better is the approximation.

Now, if you think that this article was interesting don't forget to rate it. It shows me that you care and thus I will continue write about these things.

The following two tabs change content below.
Solving ordinary differential equations (ODE)
My Twitter profileMy Facebook profileMy Google+ profileMy LinkedIn profile

Eugen Mihailescu

Founder/programmer/one-man-show at Cubique Software
Always looking to learn more about *nix world, about the fundamental concepts of math, physics, electronics. I am also passionate about programming, database and systems administration. 16+ yrs experience in software development, designing enterprise systems, IT support and troubleshooting.
Solving ordinary differential equations (ODE)
My Twitter profileMy Facebook profileMy Google+ profileMy LinkedIn profile

Latest posts by Eugen Mihailescu (see all)

Related posts:

  1. Differential equations: estimating the time of death
  2. Math can be sweet
  3. Derivative applications in real life
  4. Equation of the tangent to the curve at a point
Tagged on: ,
By Eugen Mihailescu | June 28, 2014 | Calculus, Discrete mathematics, Maths | No Comments |

Leave a Reply

Your email address will not be published. Required fields are marked *