Chapter 5 Initial-Value Problems for Ordinary Differential Equations
Updated December 9, 2003
Chapter 5 Initial-Value Problems for Ordinary Differential Equations In this chapter we look at ways to approximate solutions of differential equations. For scalar ODE, we have d y ÅÅÅÅ ÅÅÅÅÅÅ = f Hx, y, y', …, yHn-1L L, x œ I d xn n
yHx0 L = y0 , y' Hx0 L = y1 , … , yHn-1L Hx0 L = yn-1 , x0 œ I . For first order systems of ordinary differential equations we have: dY ÅÅÅÅ ÅÅÅÅ = FHx, YL dx
YHx0 L = Y0 where
y HxL y F Hx, y1 , y2 , .., yn L y jij 1 zyz jij 1 zyz jij 1 z jj y2 HxL zz jj y2 zz jj F2 Hx, y1 , y2 , .., yn L zzz jj zz jj zz jj zz YHxL = jj z , Y0 = jj zz, FHx, YL = jj zz. jj ª zzz jj ª zz jj zz ª jj zz jj zz jj zz y HxL y F Hx, y , y , .., y L k n { k n{ k n 1 2 n {
These problems are called initial-value problems (IVP).
Remark: An nth order scalar differential equation yHnL = FHx, y, y', ..., yHn-1L L can be written as a first order system of dY differential equations ÅÅÅÅ ÅÅÅÅ = GHx, YL: dx y' HxL ij yHxL yz ij yz jj zz jj zz jj y ' HxL zz jj zz y '' HxL dY j z j zz fl ÅÅÅÅ zz , GHx, YL = jj Y = jjj ÅÅÅÅ = GHx, YL zz dx z j ª ª jj zz jj zz jj zz jj zz Hn-1L Hn-1L HxL { L{ ky k FHx, y, ..., y
5.1 Elementary Theory of Initial-Value Problems
Chapter5.nb
2
We first state a theorem that guarantees that a unique solution of an IVP. dY Theorem: (Existence/Uniqueness) Consider the IVP problem: ÅÅÅÅ ÅÅÅÅ = FHx, YL , x œ I , ú YHx0 L = Y0 . Here F is a dx ∑F n+1 vector-valued function defined in a domain D of ! . Suppose F, ÅÅÅÅ ÅÅÅÅÅ , j = 1, 2, …, n, are continuous on D and ∑y j
Hx0 , Y0 L œ D . Then there exists and an open interval I such that x0 œ I on which a unique solution of the IVP exists.
Remark: It is possible to weaken the hypothesis on F. Namely, we can replace the continuity conditions on the partial derivatives of F with respect to y j by requiring that it satisfy a Lipschitz condition. That is, FHx, YL is said to satisfy a Lipschitz condition with Lipschitz constant L if » FHx, Y1 L - FHx, Y2 L » § L » Y1 - Y2 » for all Hx, Y1 L and Hx, Y2 L in each rectangle R Õ D . Here » ÿ » denotes a vector norm. In the case when n = 1 , the Existence/Uniqueness Theorem is easy to apply. Theorem: (Existence/Uniqueness). Consider the IVP dy ÅÅÅÅ d ÅxÅÅ = f Hx, yL, x œ I yHx0 L = y0
∑f If f Hx, yL and ÅÅÅÅ ÅÅÅÅ are continuous function on an open rectangle R that contains the point Hx0 , y0 L , then there exist a ∑y number h > 0 , such that there exists a unique solution of the IVP on the interval Hx0 - h, x0 + hL .
dy Example: Show that the IVP, ÅÅÅÅ ÅÅ = x y2 such that yH0L = 1 has a unique solution in a neighborhood of x = 0. dx
∑f 2 Here f Hx, yL = x y2 and ÅÅÅÅ ∑yÅÅÅ = 2 x y . These functions are continuous on ! and hence, the Existence/Uniqueness Theorem states that a unique solution will exist on some interval H-h, hL containing the origin.
Example: Consider the IVP above: dy 2ê3 ÅÅÅÅ d ÅxÅÅ = 2 y
yH0L = 0 . Does the IVP have a unique solution?
∑f -1ê3 Here f Hx, yL = 3 y 2ê3 which is continuous in all of !2 , but ÅÅÅÅ which is discontinuous in any rectangle that contains ∑ ÅyÅÅÅ = 2 y the x-axis. Since the initial condition is on the x-axis, there is no rectangle that contains the initial condition and on which f ∑f and ÅÅÅÅ ∑ ÅyÅÅÅ are continuous. Therefore, we suspect that this IVP may not have a unique solution. In fact, it doesn't!
Chapter5.nb
3
5.2 Euler's Method Consider a first order, scalar IVP: y' = f Hx, yL , yHx0 L = y0 . In general, we will not be able to find the general solution of y' = f Hx, yL and hence, it will be impossible to find the solution of the IVP. Suppose we desire to have approximations to the solution of this IVP at a set of equally spaced points: xn = x0 + n h where h is some given quantity. We have yHxL at x = x0 , namely, yHx0 L = y0 from the initial condition. Let us try to find yHx1 L. Suppose the yHxL is sufficiently smooth that we can compute its Taylor series expansion about x0 : yHxL = yHx0 L + y' Hx0 L Hx - x0 L + ÅÅÅÅ12 y'' HxL Hx - x0 L2
where x is some point between x and x0 . Let us ignore the remainder to obtain the approximation:
yHxL º yHx0 L + y' Hx0 L Hx - x0 L ï yHx1 L º yHx0 L + y ' Hx0 L Hx1 - x0 L = yHx0 L + h f Hx0 , y0 L = y1 .
The value of y1 is computable from x0 and y0 . Now let us try to find an approximation for yHx2 L. We repeat the argument: yHxL = yHx1 L + y' Hx1 L Hx - x1 L + ÅÅÅÅ12 y'' HxL Hx - x1 L2 º yHx1 L + y ' Hx1 L Hx - x1 L
ï yHx2 L º yHx1 L + y' Hx1 L Hx2 - x1 L = yHx1 L + f Hx1 , yHx1 LL Hx2 - x1 L º y1 + h f Hx1 , y1 L.
It is clear that we can continue this process and we generate approximations to yHxL at x = x1 , x2 , … . All of the approximations come from the iteration: yn = yn-1 + h f Hxn-1 , yn-1 L, n = 1, 2, … .
This is called Euler's Method. Example: Using Euler's Method, generate approximations to the solution of y' = sinHx + y2 L, yH0L = 1 for x = 0.1, 0.2, …, 1.0 . Here f Hx, yL = sinHx + y2 L, h = 0.1 , y0 = 1 , and xk = 0.1 k , k = 0, 1, 2, …, 10 .
Clear@x, y, fD; f@x_, y_D := Sin@x + y2 D; h = 1 ê 10; x@n_D := N@n * hD; y@0D = 1; y@n_ ê; n > 0D := y@nD = y@n - 1D + h * f@x@n - 1D, y@n - 1DD; Do@Print@"n = ", i, " xn = ", x@iD, " yn = ", y@iDD, 8i, 0, 10 0D := y2@nD = y2@n - 1D + h * Hx@n - 1D^ 2 + y2@n - 1D ^ 2L + Hh2 ê 2L H2 x@n - 1D ^ 2 + 2 y2@n - 1D^ 2 * Hx@n - 1D ^ 2 + y2@n - 1D^ 2LL; yrk@n_ ê; n > 0D := yrk@nD = yrk@n - 1D + Hh ê 2L * Hx@n - 1D ^ 2 + yrk@n - 1D ^ 2 + Hx@n - 1D + hL ^ 2 + Hyrk@n - 1D + h Hx@n - 1D ^ 2 + yrk@n - 1D ^ 2L ^ 2LL;
tt2 = ListPlot@Table@8x@nD, y2@nD 0D := y2@nD = y2@n - 1D + h * f@x@n - 1D, y2@n - 1DD + Hh2 ê 2L * HHD@f@xx, yyD, xxD ê. 8xx Ø x@n - 1D, yy -> f@x@n - 1D, y2@n - 1DD f@x@n - 1D, y2@n - 1DD 0D := yrk@nD = yrk@n - 1D + Hh ê 9L * H2 F1@nD + 3 F2@nD + 4 F3@nDL; tt2 = ListPlot@Table@8x@nD, y2@nD 1D := y@nD = y@n - 1D + h * f@x@n - 1D, y@n - 1DD; Do@Print@"k = ", k, " xk = ", N@x@kDD, " yk = ", N@y@kDDD, 8k, 0, 10 0 . Hence, there is no restriction on the step-size when we use this method. The phenomena were the behavior of the transient solution is much different than that of the steady-state solution is called stiffness. When a problem is stiff, then a variable step-size is probably necessary so that h can be changed to meet the conditions when the solution is changing rapidly and when the solution is slowly changing. dy Example: Use Euler's Method to approximate yH1L where ÅÅÅÅ ÅÅ = l y , yH0L = 1 with l = 100 for h = 0.1 and h = 0.01 . dx
Clear@euler, x, y, f, nD; euler@h_, k_D := Module@8 0D := ye@nD = ye@n - 1D + h * f@x@n - 1D, ye@n - 1DD; Table@8x@iD, ye@iD
Chapter 5 Initial-Value Problems for Ordinary Differential Equations In this chapter we look at ways to approximate solutions of differential equations. For scalar ODE, we have d y ÅÅÅÅ ÅÅÅÅÅÅ = f Hx, y, y', …, yHn-1L L, x œ I d xn n
yHx0 L = y0 , y' Hx0 L = y1 , … , yHn-1L Hx0 L = yn-1 , x0 œ I . For first order systems of ordinary differential equations we have: dY ÅÅÅÅ ÅÅÅÅ = FHx, YL dx
YHx0 L = Y0 where
y HxL y F Hx, y1 , y2 , .., yn L y jij 1 zyz jij 1 zyz jij 1 z jj y2 HxL zz jj y2 zz jj F2 Hx, y1 , y2 , .., yn L zzz jj zz jj zz jj zz YHxL = jj z , Y0 = jj zz, FHx, YL = jj zz. jj ª zzz jj ª zz jj zz ª jj zz jj zz jj zz y HxL y F Hx, y , y , .., y L k n { k n{ k n 1 2 n {
These problems are called initial-value problems (IVP).
Remark: An nth order scalar differential equation yHnL = FHx, y, y', ..., yHn-1L L can be written as a first order system of dY differential equations ÅÅÅÅ ÅÅÅÅ = GHx, YL: dx y' HxL ij yHxL yz ij yz jj zz jj zz jj y ' HxL zz jj zz y '' HxL dY j z j zz fl ÅÅÅÅ zz , GHx, YL = jj Y = jjj ÅÅÅÅ = GHx, YL zz dx z j ª ª jj zz jj zz jj zz jj zz Hn-1L Hn-1L HxL { L{ ky k FHx, y, ..., y
5.1 Elementary Theory of Initial-Value Problems
Chapter5.nb
2
We first state a theorem that guarantees that a unique solution of an IVP. dY Theorem: (Existence/Uniqueness) Consider the IVP problem: ÅÅÅÅ ÅÅÅÅ = FHx, YL , x œ I , ú YHx0 L = Y0 . Here F is a dx ∑F n+1 vector-valued function defined in a domain D of ! . Suppose F, ÅÅÅÅ ÅÅÅÅÅ , j = 1, 2, …, n, are continuous on D and ∑y j
Hx0 , Y0 L œ D . Then there exists and an open interval I such that x0 œ I on which a unique solution of the IVP exists.
Remark: It is possible to weaken the hypothesis on F. Namely, we can replace the continuity conditions on the partial derivatives of F with respect to y j by requiring that it satisfy a Lipschitz condition. That is, FHx, YL is said to satisfy a Lipschitz condition with Lipschitz constant L if » FHx, Y1 L - FHx, Y2 L » § L » Y1 - Y2 » for all Hx, Y1 L and Hx, Y2 L in each rectangle R Õ D . Here » ÿ » denotes a vector norm. In the case when n = 1 , the Existence/Uniqueness Theorem is easy to apply. Theorem: (Existence/Uniqueness). Consider the IVP dy ÅÅÅÅ d ÅxÅÅ = f Hx, yL, x œ I yHx0 L = y0
∑f If f Hx, yL and ÅÅÅÅ ÅÅÅÅ are continuous function on an open rectangle R that contains the point Hx0 , y0 L , then there exist a ∑y number h > 0 , such that there exists a unique solution of the IVP on the interval Hx0 - h, x0 + hL .
dy Example: Show that the IVP, ÅÅÅÅ ÅÅ = x y2 such that yH0L = 1 has a unique solution in a neighborhood of x = 0. dx
∑f 2 Here f Hx, yL = x y2 and ÅÅÅÅ ∑yÅÅÅ = 2 x y . These functions are continuous on ! and hence, the Existence/Uniqueness Theorem states that a unique solution will exist on some interval H-h, hL containing the origin.
Example: Consider the IVP above: dy 2ê3 ÅÅÅÅ d ÅxÅÅ = 2 y
yH0L = 0 . Does the IVP have a unique solution?
∑f -1ê3 Here f Hx, yL = 3 y 2ê3 which is continuous in all of !2 , but ÅÅÅÅ which is discontinuous in any rectangle that contains ∑ ÅyÅÅÅ = 2 y the x-axis. Since the initial condition is on the x-axis, there is no rectangle that contains the initial condition and on which f ∑f and ÅÅÅÅ ∑ ÅyÅÅÅ are continuous. Therefore, we suspect that this IVP may not have a unique solution. In fact, it doesn't!
Chapter5.nb
3
5.2 Euler's Method Consider a first order, scalar IVP: y' = f Hx, yL , yHx0 L = y0 . In general, we will not be able to find the general solution of y' = f Hx, yL and hence, it will be impossible to find the solution of the IVP. Suppose we desire to have approximations to the solution of this IVP at a set of equally spaced points: xn = x0 + n h where h is some given quantity. We have yHxL at x = x0 , namely, yHx0 L = y0 from the initial condition. Let us try to find yHx1 L. Suppose the yHxL is sufficiently smooth that we can compute its Taylor series expansion about x0 : yHxL = yHx0 L + y' Hx0 L Hx - x0 L + ÅÅÅÅ12 y'' HxL Hx - x0 L2
where x is some point between x and x0 . Let us ignore the remainder to obtain the approximation:
yHxL º yHx0 L + y' Hx0 L Hx - x0 L ï yHx1 L º yHx0 L + y ' Hx0 L Hx1 - x0 L = yHx0 L + h f Hx0 , y0 L = y1 .
The value of y1 is computable from x0 and y0 . Now let us try to find an approximation for yHx2 L. We repeat the argument: yHxL = yHx1 L + y' Hx1 L Hx - x1 L + ÅÅÅÅ12 y'' HxL Hx - x1 L2 º yHx1 L + y ' Hx1 L Hx - x1 L
ï yHx2 L º yHx1 L + y' Hx1 L Hx2 - x1 L = yHx1 L + f Hx1 , yHx1 LL Hx2 - x1 L º y1 + h f Hx1 , y1 L.
It is clear that we can continue this process and we generate approximations to yHxL at x = x1 , x2 , … . All of the approximations come from the iteration: yn = yn-1 + h f Hxn-1 , yn-1 L, n = 1, 2, … .
This is called Euler's Method. Example: Using Euler's Method, generate approximations to the solution of y' = sinHx + y2 L, yH0L = 1 for x = 0.1, 0.2, …, 1.0 . Here f Hx, yL = sinHx + y2 L, h = 0.1 , y0 = 1 , and xk = 0.1 k , k = 0, 1, 2, …, 10 .
Clear@x, y, fD; f@x_, y_D := Sin@x + y2 D; h = 1 ê 10; x@n_D := N@n * hD; y@0D = 1; y@n_ ê; n > 0D := y@nD = y@n - 1D + h * f@x@n - 1D, y@n - 1DD; Do@Print@"n = ", i, " xn = ", x@iD, " yn = ", y@iDD, 8i, 0, 10 0D := y2@nD = y2@n - 1D + h * Hx@n - 1D^ 2 + y2@n - 1D ^ 2L + Hh2 ê 2L H2 x@n - 1D ^ 2 + 2 y2@n - 1D^ 2 * Hx@n - 1D ^ 2 + y2@n - 1D^ 2LL; yrk@n_ ê; n > 0D := yrk@nD = yrk@n - 1D + Hh ê 2L * Hx@n - 1D ^ 2 + yrk@n - 1D ^ 2 + Hx@n - 1D + hL ^ 2 + Hyrk@n - 1D + h Hx@n - 1D ^ 2 + yrk@n - 1D ^ 2L ^ 2LL;
tt2 = ListPlot@Table@8x@nD, y2@nD 0D := y2@nD = y2@n - 1D + h * f@x@n - 1D, y2@n - 1DD + Hh2 ê 2L * HHD@f@xx, yyD, xxD ê. 8xx Ø x@n - 1D, yy -> f@x@n - 1D, y2@n - 1DD f@x@n - 1D, y2@n - 1DD 0D := yrk@nD = yrk@n - 1D + Hh ê 9L * H2 F1@nD + 3 F2@nD + 4 F3@nDL; tt2 = ListPlot@Table@8x@nD, y2@nD 1D := y@nD = y@n - 1D + h * f@x@n - 1D, y@n - 1DD; Do@Print@"k = ", k, " xk = ", N@x@kDD, " yk = ", N@y@kDDD, 8k, 0, 10 0 . Hence, there is no restriction on the step-size when we use this method. The phenomena were the behavior of the transient solution is much different than that of the steady-state solution is called stiffness. When a problem is stiff, then a variable step-size is probably necessary so that h can be changed to meet the conditions when the solution is changing rapidly and when the solution is slowly changing. dy Example: Use Euler's Method to approximate yH1L where ÅÅÅÅ ÅÅ = l y , yH0L = 1 with l = 100 for h = 0.1 and h = 0.01 . dx
Clear@euler, x, y, f, nD; euler@h_, k_D := Module@8 0D := ye@nD = ye@n - 1D + h * f@x@n - 1D, ye@n - 1DD; Table@8x@iD, ye@iD
