Resolution of first order differential equations. Method of Euler

 

Let's consider a differential equation defined by dy/dx = f(x,y) with the initial value of y(0) = y₀. Let's say there isn't any method that can't solve this equation. Not knowing how to find the solution, let's then approximate it.

Let's consider two approximate points solutions of the differential equation and draw a line through these points ( see figure above). Let's calculate the slope of that line. The slope is the tangent of the angle made by by the line and the parallel to the x-axis. Let's call this angle θ. We have:

tanθ. = yⁱᵢ₊₁ - yᵢ /xᵢ₊₁-xᵢ . It is rise/run, the formula to calculate the slope.

tanθ is also the slope of the point (xᵢ, yᵢ). Then we can write f(xᵢ, yᵢ) = yⁱᵢ₊₁ - yᵢ /xᵢ₊₁-xᵢ . We can derive the formula of Euler from this but it will not be fair to do that. What I want to show here is is a reminder of the formula of the slope which is rise/run. 

Let's consider the following figure:

We still have the same differential equation dy/dx = f(x,y) with the initial condition y(0) = y₀

Let's draw a tangent line through the point (0, y₀). Knowing x₁, let y₁ be the approximate solution of the differential equation.. Its value can be calculated using the slope of the tangent line. This slope is the slope at the point (x₀, y₀) which is f(x₀, y₀). We have two points here (x₀, y₀) and (x₁, y₁). Using the formula of the slope rise/run we have:

f(x₀, y₀) = y₁ - y₀/x₁ - x₀

(fx₀, y₀)(/x₁ - x₀) = y₁ - y₀

y₁ = y₀ + (fx₀, y₀)(/x₁ - x₀). This is the approximate value of y₁ 

Now let's draw a line through the point (x₁, y₁) and let's y₂ be the second approximation of the solution knowing x₂

The slope of the line that passes through the points (x₁, y₁) and (x₂, y₂) is :

f(x₁, y₁) = y₂ -y₁/x₂ - x₁

   f(x₁, y₁)(x₂ - x₁) = y₂ -y₁

y₂ = y₁ +  f(x₁, y₁)(x₂ - x₁). This is the approximate value of y₂.

A pattern appears here based on the values of  y₁ and y₂. If we were to calculate a third approximation of the solution which is y₃, it would be: y₃ = y₂ + f(x₂, y₂)(x₃ - x₂)'

Any approximation yₙ can be calculated by yₙ = yₙ₋₁ + f(xₙ₋₁, yₙ₋₁) (xₙ - xₙ₋₁)

Let's write  xₙ - xₙ₋₁ = h, the above equation can be written as : yₙ = yₙ₋₁ + f(xₙ₋₁, yₙ₋₁)h where h is called the step size. This is the Euler formula that allows to calculate the approximate values of a function.

Example

Let's consider a differential equation with a given solution y' + 2y = 2 - e^-4t  y(0) = 1 Use the Euler method with a step size of h = 0.1 to find approximate values of the solutions at t = 0.1 t = 0.2  t = 0.3

Solution 

Let's write the equation on its explicit form: 

y' = 2 - e^-4t - 2y

Let's apply the Euler formula:   yₙ = yₙ₋₁ + f(xₙ₋₁, yₙ₋₁)h. Since the independent variable is t, we can write the formula as yₙ = yₙ₋₁ + f(tₙ₋₁, yₙ₋₁)h

We have h = 0.1 we have to find y for t = 0.1. Since this is the first value of y, we have to call it y₁Let's substitute n by 1 in the formula:

y₁ = y₀ + f(t₀, y₀).0.1

We already have y₀. Let's calculate f(t₀, y₀) or f(0,1). In order to do that, we substitute t by 0 and y by 1 in the explicit form of the differential equation

f(0,1) = 2 - e^0 - 2 = 2 - 1 - 2 = -1

Let's substitute y₀ and f(0,1):

y₁ = 1 + (-1)(0.1) = 1 - 0.1 = 0.90

Let's find y for the value of t = 0.2: In this case we have to substitute n by 2 in the Euler formula in order to find y₂. 

y₂ = y₁ + f(t₁, y₁).0.1

Let's find f(t₁, y₁):

f(0.1, 0.9) = 2 - e^-0.4 - 2(0.9) = -0.47

y₂ = 0.9 + (-0.47)(0.1) = 0.85.

Let's find another approximate value of y for t = 0.3

y₃ = y₂ + f(t₂,y₂).0.1

Let's calculate f(t₂,y₂):

f(0.2, 0.85) = 2 - e^-0.8 -2.0.85 = 2 - 1/e^0.8 - 1.70 = 2 - 0.449 - 1.70 = -0.149

y₃ = 0.85 + (-0.149)(0.1)

y₃ = 0.8351

Where the calculations of the approximate values contain decimals it's better to find the approximations with more than 2 or 3 decimals.

Interested in reviewing Calculus I check this site: Center for Integral Development 

Interval of validity of a non-linear differential equation

 In the following problem we are going to determine the interval of validity of a non-linear differential equation.

Problem. Determine the interval of validity of the initial value problem and give its dependence on the value of y₀.

y' = y²  y(0) = y₀

Let's see if f(t, y) = y² and δf/δy are continuous. It's obvious that f(t, y) =  y² is continuous. 

Let's calculate δf/δy:

δf/δy = 2y. It's continuous. Then there is a unique solution to the IVP in an interval that contains t₀ = 0

We have to determine this solution in order to find the interval of validity. 

Let's solve the differential equation

dy/dt = :y² 

dy = y² dt

dy/ y² = dt

y^-2 dy = dt

∫y^-2 dy = ∫dt

y^-2+1/-2 +1 = t + c

y^-1/-1 = t

-1/y = t + c

y = -1/t + c

Let's determine the value of c. Let's substitute t by 0 and y by  y₀ as it's given in the initial condition

y₀ = -1/c c = -1/y₀

Let's substitute c in y we have:

y = -y₀ /y₀ t -1

The interval of validity of the IVP is the interval where y is defined. 

Let's determine this interval :of validity.

The function y isn't defined if the denominator is equal to 0.

If y₀ t -1 = 0, t = 1/y₀ . Therefore y isn't defined for that value of t.

Y is defined in the following intervals: -∞<t<1/y₀  and 1//y₀<t<+∞

Let's determine the interval of validity for different values of y₀:

If y₀ = 0, then the interval of validity is -∞<t<+∞

If y₀<0  1/y₀<0 t = 0 doesn't belong to  -∞<t<1/y₀ but belongs to 1//y₀<t<+∞, which is an interval of validity

If y₀>0 then 1/y₀>0  -∞<t<1/y₀

The problem that we just solved demonstrates the difference between the interval of validity of linear differential equations and that of non-linear differential equations. The interval of validity of linear differential equations doesn't depend the value of y₀. The interval of validity of non-linear differential equations depends on the value of  y₀.

Intervals of validity (continued)

Theorem. Let's consider the Initial Value Problem (IVP) y' = f(t,y)  and y(t₀) = y₀. If f(t,y) and ẟf/ẟy are continuous in some rectangle α<t<β and  δ<y<ϒ containing the point (t₀,y₀) then there is a unique solution to the IVP in the interval t₀ - h<t< t₀ + h that is contained in α<t<β. 

This theorem doesn't allow to find the interval of validity. If its conditions are met, we know that the unique solution exists. This unique solution will be needed to determine the interval of validity.

In non-linear differential equations the value of y₀ can affect the interval of validity.

Example. Solve the initial value problem y' = y^1/3 y(0) = 0

Solution

Let's see  if the conditions are met

f(t,y) = y^1/3 is continuous on any interval.

Let's see if the derivative is continuous. The derivative of f is 

δf/δy = (1/3)y^1/3-1

         = (1/3)y^-2/3

         =  (1/3) 1/y^2/3

         =  1/3y^2/3

The function  δf/δy isn't continuous for y = 0. The theorem states that the function must be continuous in an interval that contains  y₀. Since the function isn't continuous for y = 0, it's not continuous in any interval that contains  y₀.= 0. The conditions aren't met. There isn't a unique solution to the differential equation in an interval that contains y = 0. Let's solve this equation n order to find its solutions.

dy/dt = y^1/3

dy =   y^1/3dt

Let's separate the variables;

dy/ y^1/3 = dt

 1/y^1/3dy = dt

y^-1/3dy = dt

Let's integrate both sides:

∫ y^-1/3dy = ∫dt

y^-1/3 + 1/-1/3 + 1 = t + c

y^2/3 /2/3 = t + c

3/2y^2/3 = t + c

For t = 0 y = 0. Let's substitute t and y to find c:

0 = 0 + c

c = 0. 

Let's substitute c in the last equation:

 3/2y^2/3 = t

 y^2/3 = t/3/2

   y^2/3  = 2/3t

Let's raise both sides to the third power:

 (y^2/3)³ = (2/3t)³

y² =  (2/3t)³

y = ∓⎷ (2/3t)³

y =  ∓ (2/3t)^3/2 

These solutions satisfy the initial condition t = 0 y = 0. y(t) = 0 is a solution since it verifies the initial condition t = 0 y = 0.

If you are interested in reviewing Calculus I or know someone interested, please refer to Center for Integral Development 

 

 

Intervals of validity in differential equations

 It seems normal to solve a differential equation in order to find its interval of validity but it's not always the case. The following theorem will allow us to find the interval of validity without solving the differential equation. An interval of validity guarantees the existence of a solution and its uniqueness.

Here is the theorem:

Let's consider the following differential equation y' + p(t)y = g(t) with the initial value y(t₀) = y₀. If p(t) and g(t) are continuous on the interval α<t<β and the interval contains t₀, then there is a unique solution of the differential equation on that interval. If the interval is the largest one in which p(t) and g(t) are continuous, then this interval is the interval of validity of the equation.

Example. Determine the interval of validity of the initial value problem without solving it: 

(t² - 9)y' + 2y =ln❘20 - 4t❘  y(4) = -3

Let's write the equation under its explicit form i.e y' + p(t)y = g(t). In order to do so, we divide the equation by the coefficient of y' which is t² - 9:

(t² - 9)y'/(t² - 9) + 2y/(t² - 9) = ln❘20 - 4t❘/t² - 9

y' + 2y/(t² - 9) =  ln❘20 - 4t❘/t² - 9

Let's study the continuity of the functions p(t) = 2/ t² - 9 and g(t) = ln❘20 - 4t❘/t² - 9.

p(t) being a rational function is not defined for the value of t that make its denominator equal to zero.

t² - 9 = 0 for t = -3 or t = 3. Therefore p(t) is continuous in the intervals ]-∞  -3[ ]-3 3[ ]3  +∞[. 

The function g(t) is continuous if  its numerator is defined and its denominator are not equal to zero. The numerator is undefined if 20 - 4t = 0 because in this case we would have ln0. The numerator is defined in the intervals ]-∞ 5[ 5[  +∞[. The denominator is defined in the intervals ]-∞ --3[]-3 3[ ]3  +∞[. The function g(t) is continuous in the union of the intervals in which both the numerator and the denominator are defined. The function g(t) is continuous in ]-∞ -3[]-3 3[ ]3 5[ ]5 +∞[. These are also intervals in which both p(t) and g(t) are continuous. The interval of validity is the interval that has the initial value t = 4. In this case this is the interval ]3 5[. 

Interested in reviewing Differential and Integral Calculus check Center for Integral Development

Resolution of Differential equations having the form y' = G(ax + by)

This is another form of differential equation that can be solved using the substitution method. Remember that the substitution method was used in the case of homogeneous differential equations. Homogeneous equations have the form y' = F(y/x).

In order to solve the differential equation  y' = G(ax + by), we use the substitution v = ax + by. Then we have y' = G(v)

Let's derive the equation   v = ax + by

v' = a + by' 

by' = v' - a

y' =  (v' - a)/b

Let's equal both expressions of y'

 (v' - a)/b = G(v)

v' - a = bG(v)

v' = a + bG(v)

dv/dx = a + bG(v)

1/dx = a + bG(v)/dv

dx = [1/a + bG(v)]dv

Let's integrate both sides:

∫dx = 1/a + bG(v∫dv

The resolution of this equation allows to find v. Once we have v we can substitute it in v = ax + by that allows to find y.

Example. Solve the differential equation: y' - (4x - y + 1)² = 0

Let's write v = 4x - y

Then the differential equation becomes: y' - (v + 1)² = 0 (1)

Let's derive v in order to obtain y'. Once we find y' we will have a differential equation in v. We will have to solve this equation for v. Once we find v we can substitute to find y.

v' = 4 - y' then y' = 4 - v'.

Let's substitute y' in (1):

4 - v' - (v + 1)² = 0

4 - (v + 1)² = v'

dv/dx = 4 - (v + 1)²

1/dx = 4 - (v + 1)²/dv

dx = dv/4 - (v + 1)²

Let's integrate both sides:

∫dx = ∫dv/4 - (v + 1)² (2)

Let's calculate the second side. In order to do this, let's transform 1/4 - (v + 1)²:

1/4 - (v + 1)² = 1/-v² -2v + 3

                      = -1/v² + 2v -3

                      = - 1/(v + 3) (v - 1)

   Then the equation (2) becomes:

∫dx = - ∫1/(v + 3) (v - 1) dv (3)

Let's develop 1/(v + 3) (v - 1) in partial fractions:

1/(v + 3) (v - 1) = A/v + 3 + B/v-1

By developing this expression we find A = -1/4 and B = 1/4

By substituting we find :

1/(v + 3) (v - 1) = -1/4( v + 3)  + 1/4(v - 1)

The equation (3) becomes:

 ∫dx  = ∫[ 1/4( v + 3)  - 1/4(v - 1)]dv

 ∫dx = 1/4∫(1/v + 3 -1/v - 1)

x + c = 1/4 ln( v + 3/v - 1) (we don't consider a constant of integration for this integral)

4x + 4c = ln( v + 3/v - 1)

v + 3/v - 1 = e^4x + 4c

                 = e^4x.e^4c

                = (e^4x.)(c) ( we consider 4c to be equivalent to another constant c which we can call also c')

v + 3 = ce^4xv - ce^4x

By doing all the calculations, we find v = c^4x + 3/ce^4x - 1

By substituting v in the expression v = 4x - y we find y = 4x - 3 + ce^4x/-1 + ce^4x

Practice. Solve the differential equation y' = e^9y-x

Interesting in reviewing the Differential and Integral Calculus check out this site Center For Integral Development 

Homogeneous Differential Equations

 An homogeneous differential equation is one that can be written as y' = F(y/x). The equation is not written in this form most of the times. You have to transform it in order to write it under this form. This type of equation is solved by substitution. That's another type of equation that's solved by substitution. Remember that in the case of a Bernoulli differential equation we use the substitution v = y¹⁻ⁿ. In the case of an homogenous differential equation we are going to see what type of substitution we can use.

Let's bring this homogeneous differential equation to us in order to solve it. Here it comes as y' = F(y/x). It seems natural that we can substitute y/x by something. Let's call it v. We can write v =  y/x. Then we have: y = vx.

Let's derive both sides of the equation:

y' = v'x + v

Remember also that the differential equation is y' = F(y/x) or y' = F(v)

 Equaling both expressions of y' we have:

v'x + v = F(v)

We have a differential equation with a new dependent variable v. This equation can be solved by the method of separation of variables. Let's continue;

v'x = F(v) - v

(dv/dx)x = F(v) - v

x/dx = F(v) - v/dv

dx/x  = dv/F(v) - v

Now that the variables are separated we can integrate both sides of the equation in order to find v. Once we find v we can use the substitution to find y.

Example. Solve the differential equation xyy' + 4x² + y² = 0

Let's bring the expression y/x in the equation by dividing by x. Let's divide both sides by x²:

 xyy'/x² + 4x²/x² +  y²/x²  = 0

yy'/x + 4 + (y/x)² = 0

y/xy' + 4 + (y/x)² = 0

Let's substitute y/x by v:

vy' + 4 + v² = 0  (1)

Let's find y' in order to have an equation with the new variable v:

y/x = v 

y = vx

y' = v'x + v

Let's substitute y' in equation (1)

v(v'x + v) + 4 + v² = 0

vv'x + v² + 4 +  v² = 0

vv'x + 2v² + 4 = 0

vv'x  = -2v² - 4 

Let's divide both sides by v in order to separate the variables:

v'x  = -2v² - 4/v

(dv/dx)x = -2v² - 4/v

x/dx = -2v² - 4/vdv

dx/x = -vdv/ 2v² + 4

Now that the variable are separated let's integrate both sides in order to find v;

∫dx/x = ∫(-v/ 2v² + 4)dv = -∫(v/ 2v² + 4)dv (2)

Let's solve the integral in the second side by doing the substitution u =  2v² +4

du = 4vdv dv = du/4v

The integral at the second side becomes:

∫(v/ 2v² +4)dv = ∫(v/u)(du/4v) = ∫du/4u = 1/4∫du/u = 1/4lnu = 1/4ln (2v² +4) + k

Let's substitute the value of the integral we just calculated in (2)

∫dx/x = -1/4ln (2v² + 4) + k

lnx + K' = -1/4ln (2v² +4) + k

1/4ln (2v² + 4) = -lnx - k' + k

Let's call c the expression-k'+k 

1/4ln (2v² + 4) = -lnx + c

1/4ln (2v² + 4) = -lnx + c

ln(2v² + 4)^1/4  =  -lnx + c

(2v² + 4)^1/4 = e^(-lnx + c)

                     = (e^-lnx)(e^c)

                     = (e^lnx^-1)(d) We call e^c d

  e^lnx^-1 being equal to x^-1, we have:

(2v² + 4)^1/4 = (x^-1)(d) 

Let's raise both sides to the power of 4:

2v² + 4 = x⁻⁴d⁴

Let's call d⁴ g:

2v² + 4 = x⁻⁴g

            = g/x⁴

v²  = g-4x⁴/2x⁴

v = ∓(⎷2/2) (⺁g-4x⁴)/x²

By substituting v in the expression v = y/x, we obtain:

y = ∓⎷2/2⺁g-4x⁴/x

Practice. Solve the differential equation for the initial condition y(2) = -7 and x>0

Bernoulli Differential Equations

  Let's consider a differential equation in the form y' + p(x)y = q(x)yⁿ where p(x) and q(x) are continuous and n a real number. Such an equation is called Bernoulli equation. If n = 0 or n = 1, the equation becomes easy to solve as an ordinary differential equation.

We are going to use the substitution method to solve this equation. It appears that the term yⁿ is a problem. Let's get rid of it by multiplying both sides by its inverse.

y⁻ⁿy' +y⁻ⁿp(x)y =  q(x)yⁿy⁻ⁿ

y⁻ⁿy' + p(x)y¹⁻ⁿ = q(x)

Now let's do the substitution that will lead us to an equation easier to solve. Let's write v = y¹⁻ⁿ 

Let's differentiate both sides. v is a function of y and y is a function of x. We use the derivative formula of uⁿ = nuⁿ⁻¹u' to find the derivative of  y¹⁻ⁿ 

v' = (1-n) y¹⁻ⁿ⁻¹y'

y' = v'/(1-n) y¹⁻ⁿ⁻¹ = v'/(1-n)y⁻ⁿ

Now if we substitute y' and v in the equation we will have an equation in v easier to solve. Let's do that.

y⁻ⁿv'/(1-n)y⁻ⁿ + p(x)v = q(x)

Let's simplify the first expression by y⁻ⁿ:

v'/(1-n) + p(x)v = q(x)

This is a linear differential equation in v that can be solved by finding v. Once we find v we can substitute it and find y. Let's do a numerical equation.

Example

Let's solve the equation: y' + (4/x)y = x³y² 

Let's multiply both sides by y⁻² to get rid of y² :  y⁻²y' + (4/x)yy⁻² = x³y² y⁻² (1)

y⁻²y' + (4/x)y⁻¹ = x³ (2)

Let's write v = y⁻¹

Let's derive both sides:

v' = (-1)y⁻¹⁻¹y'

v' = -y⁻²y'

Let's find y': 

y' = v'/-y⁻²

y' = - v'y²

Let's substitute v and y' in the equation (2):

(y⁻²)( - v'y²) + (4/x)v = x³

- v' +  (4/x)v = x³

v' -(4/x)v = -x³ (3)

To solve this linear differential equation in v we have to find μ:

μ = e^∫p(x)dx

p(x) is the coefficient of v: p(x) = -4/x

Let's calculate ∫p(x)dx

∫p(x)dx = ∫-4/xdx = -4lnx = lnx⁻⁴

Let's substitute ∫p(x)dx  in the formula of μ:

μ = e^ lnx⁻⁴ = x⁻⁴

Let's multiply both sides of equation (3) by x⁻⁴:

v'x⁻⁴ - x⁻⁴.(4/x)v = (x⁻⁴)( -x³)

v'x⁻⁴ - 4x⁻⁵v = -x⁻¹

The left side is the derivative of (x⁻⁴v):

(x⁻⁴v)' = -x⁻ ¹

Let's use the notation d/dx:

d(x⁻⁴v)/dx = -x⁻¹

d(x⁻⁴v) = -x⁻¹dx

Let's integrate both sides:

∫d(x⁻⁴v) = ∫-x⁻¹dx

x⁻⁴v = -∫1/xdx = -lnx + k

v = (-lnx + k)/x⁻⁴

v = x⁴(-lnx + k)

Let's substitute v in v = y⁻¹:

 x⁴(-lnx + k) = y⁻¹:

y = 1/x⁴(-lnx + k)

This value exists only if x is different of 0.

Interested in reviewing Differential and Integral Calculus visit Center For Integral Development