Download Introduction to ODEs - Mathematical Modelling - Notes | MATH 456 and more Study notes Mathematics in PDF only on Docsity!
Chapter 2
Introduction to ODEs
In the following short chapter we will summarize some easy and efficient methods for solving simple first order Ordinary Differential Equations (ODEs). In order to use the appropriate method to obtain the solution of a given ODE we must be able to classify the type of ODE we have at hand.
2.1 Classification
In general we classify differential equations as either partial or ordinary. For now we examine only ordinary differential equations. In general however most models are usually comprised of several (a system) partial or ordinary differential equations. Another reason that ordinary differential equations (ODEs) are important is that techniques of solution for partial differential equations (PDEs) can stem from those in ordinary differential equations. So learning how to solve ODEs will be useful by itself or in conjuction to solving PDEs. To choose the appropriate method of solution for a given ODE you must find out what type of ODE you have. ODEs are categorized, among other things, as: linear/non-linear, ordinary/partial. Also we usually refer to the “order” of the ODE as well as to whether it is “homogeneous” or not. The following examples should be helpful in understanding how this classification system works:
o dy dx = − 3 x + 5 1st order linear non-homogeneous ODE o 3 dy dx = 9y^2 1st order non-linear homogeneous ODE o 3 dy dx + y = 0 1st order linear homogeneous ODE o 3 d
(^2) y dx^2 + 9x
(^3) y = 9 2nd order linear non-homogeneous ODE
o 7
(dy dx
- 4y = 0 1st order non-linear homogeneous ODE
Try to figure out based on the examples above how naming of ODEs works. You probably understand from the above how we determine the order of an ODE - based on the order of the derivative in the equation. How about linearity? Do you see it? Linearity refers to the power of y. For instance anything with y^2 and higher is non-linear. Similarly if a derivative was raised to a power, such as ( dy dx )^2 , then it is also non-linear. Now determining whether an ODE is homogeneous or not is a different story. The idea, as you may have noticed, is to bring all terms which include y’s on one side of the equation. If by doing so the other side of the equation is 0 then you have a homogeneous ODE. Otherwise you have a non-homogeneous ODE.
2.2 Solving first order homogeneous ODEs
Now that you have learned some of the basics about the names of these equations you will put them to use in order to decide which method to choose in order to obtain their solutions. If you correctly identified the ODE then you will be able to solve it! It is not the point of our class to review all methods of solution for a given ODE. Instead we emphasize on the ones we will probably put to use during our study of mathematical models. In that respect we start by studying methods of solution for 1st order homogeneous ODEs.
2.2.1 Non-linear first order homogeneous ODEs
It may sound surpricing but some of the easiest first order ODEs to solve are the nonlinear (homoge- neous) ones. Let us examine for example the following non-linear 1st order homogeneous differential equation,
xy ′ = y^2 , or otherwise also written as: x
dy dx
= y^2
with the following initial condition y(1) = −1. We will present the solution of this equation based on the separation of variables method. We outline this method of solution for such an ODE below:
- Separate variables: dy y^2
dx x
- Integrate both sides: (^) ∫ dy y^2
dx =
dx x
dx
which gives, −
y
= ln(x) + C (2.1)
do not forget the constant of integration C! To evaluate C we use the initial condition, y(1) = 5. If, on the other hand, no initial condition is provided then C will be part of the final solution for this equation. Note that the initial condition y(1) = 5 implies that for x = 1 then y = 5. Substituting these values to equation (2.1) above we obtain,
C = −
− ln(1) = −. 2
Thus equation (2.1) becomes −
y
= ln(x) −. 2
That’s it! We solved this differential equation! How do we know that we solved it? Easy! Do you see any derivatives left? If not then we are done. The equation is solved! In this case we can even go further and solve it for y itself (something that is not always necessary or sometimes even possible), y = −
ln(x) − (^15)
Ok, so now we can take care of the non-linear homogeneous type first order ODEs. Let’s now see how we can solve the linear ones.
2.3 Numerical solutions to ODEs
In this section we will learn how to solve a number of different types of ODEs using a computer and applying some well known and very effective numerical algorithms. We start our exposition with a simple but robust method.
2.3.1 Euler’s method
This method is simple to learn and will essentially allow us to solve almost any type of ODE which we will come across. Naturally we implement this method in the computer but in fact for simple problems you could possibly obtain the solution with a calculator (or even by hand). The starting point is to write your ODE in the following form,
y
′ = f (x, y)
where f (x, y) corresponds to any other term in your ODE. For instance looking back at the equation xy ′ = y^2 which we solved earlier we would rewrite it as,
y
′
y^2 x
In this case therefore f (x, y) = y
2 x. Once this is established the method iterates in the computer with the following formula, Euler’s method: yn+1 = yn + hf (xn, yn)
where here h denotes the step size in x. The iteration counter n starts at 0 and will go on until we produce the solution required by the problem. For instance recall that the initial condition for our ODE was provided to be y(1) = 5. In this case for n = 0 we have that x 0 = 1 and y 0 = 5 with which to start the iteration of our method. The only thing which may be left up to you (unless specifically given by the problem) is how big the step size h should be chosen to be. Euler’s method is not exact! In fact it almost always includes errors and the solution predicted is only an approximation to the true solution. The bigger the step size h the bigger the errors to our solution. So if you would like to have a very accurate solution then you must take a very small step size h and therefore you must iterate several times. Let us assume that you would like to obtain the solution y when x = 10 for the following ODE
xy ′ = y^2 , with initial condition y(1) = − 1
Remember that we rewrite the ODE first. So here f (x, y) = y^2 /x. Since the step size h is not specified you are free to choose. So if you take h = 1 then you would have to iterate Euler’s method 9 times since you are starting at x = 1 and need to reach x = 10 with this h. The following table shows the result of the computer output and how each y is produced after each iteration of the method,
n 0 1 2 3 4 5 6 7 8 9 x 1 2 3 4 5 6 7 8 9 10 y − 1 −. 5 −. 42 −. 37 −. 34 −. 32 −. 31 −. 29 −. 28 −. 2801 f (x, y). 5. 08. 04. 03. 02. 01. 01. 01. 008. 007
So the solution at x = 10 is y = −.2801 using a step size of h = 1. Is this correct? In fact we can check this since we have found earlier the exact solution to be
y = −
ln(x) + 1
Let us compare, side by side, the approximate solution produced by Euler’s method above with this exact solution in the table below:
x 1 2 3 4 5 6 7 8 9 10 yapp. − 1 −. 50 −. 41 −. 37 −. 34 −. 32 −. 3104 −. 2983 −. 2884 −. 2801 yex − 1 −. 59 −. 47 −. 41 −. 38 −. 35 −. 3395 −. 3247 −. 3128 −. 3028
This gives that y = −.2801 for x = 10 while the equivalent exact solution is y = .3028. The absolute error between the two being .0227 and the relative error .075. We can not help but wonder how much better we could possibly do if we decrease the step size h = .5. In that case of course it would take us 18 iterations starting at x = 1 to reach x = 10. We display some of these results in the table below:
n 0 1 2 3 4... 15 16 17 18 x 1 1. 5 2 2. 5 3... 8. 5 9 9. 5 10 yapp. − 1 −. 66 −. 55 −. 49 −. 45... −. 3091 −. 3037 −. 2989 −. 2944 yex − 1 −. 71 −. 59 −. 52 −. 47... −. 3185 −. 3128 −. 3076 −. 3028
This time the absolute error between the exact and approximate solutions for x = 10 is .0084 while the relative error is just .000003. Huge improvement for a little bit extra work for the computer.
2.3.2 Runge-Kutta Method
We now present an even more accurate method which although it is slightly more complicated to program in the computer can also solve an ODE of the type,
y ′ = f (x, y)
Usually the Runge-Kutta method is more accurate than the Euler method for the same value of the step size h. The Runge-Kutta method goes as follows:
yn+1 = yn +
(k 1 + 2k 2 + 2k 3 + k 4 ) where k 1 = hf (xn, yn)
k 2 = hf (xn +
h, yn +
k 1 )
k 3 = hf (xn +
h, yn +
k 2 )
k 4 = hf (xn + h, yn + k 3 )
Just so that we can see the differences between the Euler and Runge-Kutta methods we solve the same exact ODE as before,
xy
′ = y^2 , with initial condition y(1) = − 1
In this format it is not that hard to see that we can write this system into a matrix as follows,
U
′ = AU where A =
In this format it is now clear that our original single 3rd order ODE (2.2) has been transformed to a reduced 1st order system of three ODEs. In general this method can take a single nth order ODE and transform it into a system of n first order ODEs. Once you have such a system you can apply Euler or Runge-Kutta method to solve it. The following simple example should be illustrative of this procedure. Example: Solve the following ODE using Euler’s method. Obtain the solution for x = 1.
y
′′ − y = x with initial conditions: y(0) = 0, y
′ (0) = 1 (2.4)
We start by reducing the 2nd order ODE into a system of 2 first order ODEs. To do this we follow the procedure outlined previously. We first define 2 new variables, u 1 and u 2 as follows,
u 1 = y u 2 = y ′
Now we differentiate both sides, u
′ 1 =^ y
′ u
′ 2 =^ y
′′
and eliminate y ′′ by solving (2.4). Thus, u ′ 2 =^ y
′′ = x + y and replace y here with u 1. Thus,
u ′ 1 =^ y
′ = u 2 = 0u 1 + 1u 2 u
′ 2 =^ x^ +^ y^ =^ x^ +^ u^1 = 1u^1 + 0u^2 +^ x
To make the notation clearer we rename u 1 = U and u 2 = V. Thus we have the following two equations, U ′ = F (x, U, V ) where F (x, U, V ) = V V ′ = G(x, U, V ) where G(x, U, V ) = x + U
We are now ready to apply Euler’s method for this system. Before we start though we should also translate the initial conditions to correspond to U and V. Starting with y(0) = 0 and since y = u 1 = U then we have that U(0) = 0. Similarly the other initial condition y ′ (0) = 1 becomes V (0) = 1 since y ′ = u 2 = V. Euler’s method for our system is equivalent to the following:
Un+1 = Un + hF (xn, Un, Vn) Vn+1 = Vn + hG(xn, Un, Vn)
The following table provides the solution assuming a step size of h = .2,
n 0 1 2 3 4 5 x 0. 2. 4. 6. 8 1 U 0. 2. 4. 616. 864 1. 1606 V 1 1 1. 08 1. 24 1. 48 1. 816
Thus when x = 1 we obtain the approximate solution to be U = u 1 = 1.16 and V = u 2 = 1.82. This solves our problem! Well not really. Remember that the original question was to solve for y not U, V, u 1 or even for u 2. However we know that y = u 1. Thus y ≈ 1 .16 when x = 1! Just to check our solution if we reduce the step size to h = .1 and iterate 10 times we obtain the solution to be y ≈ 1 .24 when x = 1. If the step size becomes h = .01 and we iterate 100 times we obtain y ≈ 1 .3388. Last if we reduce the step size to h = .001 and iterate 1000 times we obtain that y ≈ 1 .3492. In other words we can safely believe that y ≈ 1 .34.
