r""" This file contains a few routines which are useful to teaching a lower level ODE course or a more advanced level calculus I course. The goal is to find an approximate solution to the problem \begin{equation} \label{eqn:*} y'=f(x,y),\ \ \ y(a)=c, \end{equation} where $f(x,y)$ is some given function. We shall try to approximate the value of the solution at $x=b$, where $b>a$ is given. The basic idea can also be explained ``algebraically''. Recall from the definition of the derivative in calculus 1 that \[ y'(x)\cong \frac{y(x+h)-y(x)}{h}, \] $h>0$ is a given and small. This an the DE together give $f(x,y(x))\cong \frac{y(x+h)-y(x)}{h}$. Now solve for $y(x+h)$: \[ y(x+h)\cong y(x)+h\cdot f(x,y(x)). \] If we call $h\cdot f(x,y(x))$ the ``correction term'' (for lack of anything better), call $y(x)$ the ``old value of $y$'', and call $y(x+h)$ the ``new value of $y$'', then this approximation can be re-expressed \[ y_{new}=y_{old}+h\cdot f(x,y_{old}). \] {\bf Tabular idea}: Let $n>0$ be an integer, which we call the {\bf step size}. This is related to the increment by \[ h=\frac{b-a}{n}. \] This can be expressed simplest using a table. \begin{center} \begin{tabular}{c|c|c} $x$ & $y$ & $hf(x,y)$ \\ \hline $a$ & $c$ & $hf(a,c)$ \\ $a+h$ & $c+hf(a,c)$ &\vdots \\ $a+2h$ & \vdots & \\ \vdots & & \\ $b$ & ??? & xxx \\ \end{tabular} \end{center} The goal is to fill out all the blanks of the table but the xxx entry and find the ??? entry, which is the {\bf Euler's method approximation for $y(b)$}. \vskip .2in \begin{center} Improved Euler's method \end{center} {\bf Geometric idea}: The basic idea can be easily expressed in geometric terms. As in Euler's method, we know the solution must go through the point $(a,c)$ and we know its slope there is $m=f(a,c)$. If we went out one step using the tangent line approximation to the solution curve, the approximate slope to the tangent line at $x=a+h, y=c+h\cdot f(a,c)$ would be $m'=f(a+h,c+h\cdot f(a,c))$. The idea is that instead of using $m=f(a,c)$ as the slope of the line to get our first approximation, use $\frac{m+m'}{2}$. The ``improved'' tangent-line approximation at $(a,c)$ is: \[ y(a+h)\cong c+h\cdot \frac{m+m'}{2} =c+h\cdot \frac{f(a,c)+f(a+h,c+h\cdot f(a,c))}{2}. \] (This turns out to be a better apprpximation than the tangent-line approximation $y(a+h)\cong c+h\cdot f(a,c)$ used in Euler's method.) Now we know the solution passes through a point which is ``nearly'' equal to $(a+h,c+h\cdot \frac{m+m'}{2})$. We now repeat this tangent-line approximation with $(a,c)$ replaced by $(a+h,c+h\cdot f(a,c)$. Keep repeating this number-crunching at $x=a$, $x=a+h$, $x=a+2h$, ..., until you get to $x=b$. {\bf Tabular idea}: The integer step size $n>0$ is related to the increment by \[ h=\frac{b-a}{n}, \] as before. The improved Euler method can be expressed simplest using a table. \begin{center} \begin{tabular}{c|c|c} $x$ & $y$ & $h\frac{m+m'}{2}=h\frac{f(x,y)+f(x+h,y+h\cdot f(x,y))}{2}$ \\ \hline $a$ & $c$ & $h\frac{f(a,c)+f(a+h,c+h\cdot f(a,c))}{2}$ \\ $a+h$ & $c+h\frac{f(a,c)+f(a+h,c+h\cdot f(a,c))}{2}$ &\vdots \\ $a+2h$ & \vdots & \\ \vdots & & \\ $b$ & ??? & xxx \\ \end{tabular} \end{center} The goal is to fill out all the blanks of the table but the xxx entry and find the ??? entry, which is the {\bf improved Euler's method approximation for $y(b)$}. The idea for systems of ODEs is similar. This is implemented below as well. AUTHOR: David Joyner (1-2006) bug fixes -- 4-30-2006 """ #***************************************************************************** # Copyright (C) 2006 David Joyner # 2006 William Stein # # Distributed under the terms of the GNU General Public License (GPL) # # http://www.gnu.org/licenses/ #***************************************************************************** def improved_eulers_method(f,x0,y0,h,x1): """ This implements Improved Euler's method for finding numerically the solution of the 1st order ODE y' = f(x,y), y(a)=c. The "x" coloumn of the table increments from x0 to x1 by h (so (x1-x0)/h must be an integer). In the "y" column, the new y-value equals the old y-value plus the corresponding entry in the last column. EXAMPLES: sage: RR = RealField(sci_not=0, prec=4, rnd='RNDU') sage: x,y = PolynomialRing(RR,2).gens() sage: improved_eulers_method(5*x+y-5,0,1,1/2,1) x y (h/2)*(f(x,y)+f(x+h,y+h*f(x,y)) 0 1 -1.87 1/2 -0.875 -1.37 1 -2.25 -0.687 sage: x,y=PolynomialRing(QQ,2).gens() sage: improved_eulers_method(5*x+y-5,0,1,1/2,1) x y (h/2)*(f(x,y)+f(x+h,y+h*f(x,y)) 0 1 -15/8 1/2 -7/8 -95/64 1 -151/64 -435/512 AUTHOR: David Joyner """ print "%10s %30s %50s"%("x","y","(h/2)*(f(x,y)+f(x+h,y+h*f(x,y))") n=((RealField(max(6,RR.precision()))('1.0'))*(x1-x0)/h).integer_part() x00=x0; y00=y0 for i in range(n+ZZ(1)): s1 = f(x00,y00) s2= f(x00+h,y00+h*f(x00,y00)) corr = h*(ZZ(1)/ZZ(2))*(s1+s2) print "%10r %30r %30r"%(x00,y00,corr) y00 = y00+h*(ZZ(1)/ZZ(2))*(f(x00,y00)+f(x00+h,y00+h*f(x00,y00))) x00=x00+h def eulers_method(f,x0,y0,h,x1): """ This implements Euler's method for finding numerically the solution of the 1st order ODE y' = f(x,y), y(a)=c. The "x" column of the table increments from x0 to x1 by h (so (x1-x0)/h must be an integer). In the "y" column, the new y-value equals the old y-value plus the corresponding entry in the last column. EXAMPLES: sage: x,y=PolynomialRing(QQ,2).gens() sage: eulers_method(5*x+y-5,0,1,1/2,1) x y h*f(x,y) 0 1 -2 1/2 -1 -7/4 1 -11/4 -11/8 sage: RR = RealField(sci_not=0, prec=4, rnd='RNDU') sage: x,y=PolynomialRing(RR,2).gens() sage: eulers_method(5*x+y-5,0,1,1/2,1) x y h*f(x,y) 0 1 -2.00 1/2 -1.00 -1.75 1 -2.75 -1.37 AUTHOR: David Joyner """ print "%10s %20s %25s"%("x","y","h*f(x,y)") n=int((RealField(max(6,RR.precision()))('1.0'))*(x1-x0)/h) x00=x0; y00=y0 for i in range(n+ZZ(1)): print "%10r %20r %20r"%(x00,y00,h*f(x00,y00)) y00 = y00+h*f(x00,y00) x00=x00+h def eulers_method_2x2(f,g, t0, x0, y0, h, t1): """ This implements Euler's method for finding numerically the solution of the 1st order system of two ODEs x' = f(t, x, y), x(t0)=x0. y' = g(t, x, y), y(t0)=y0. The "t" coloumn of the table increments from t0 to t1 by h (so (t1-t0)/h must be an integer). In the "x" column, the new x-value equals the old x-value plus the corresponding entry in the next (third) column. In the "y" column, the new y-value equals the old y-value plus the corresponding entry in the next (last) column. EXAMPLES: sage: t, x, y = PolynomialRing(QQ,3).gens() sage: f = x+y+t; g = x-y sage: eulers_method_2x2(f,g, 0, 0, 0, 1/3, 1) t x h*f(t,x,y) y h*g(t,x,y) 0 0 0 0 0 1/3 0 1/9 0 0 2/3 1/9 7/27 0 1/27 1 10/27 38/81 1/27 1/9 sage: RR = RealField(sci_not=0, prec=4, rnd='RNDU') sage: t,x,y=PolynomialRing(RR,3).gens() sage: f = x+y+t; g = x-y sage: eulers_method_2x2(f,g, 0, 0, 0, 1/3, 1) t x h*f(t,x,y) y h*g(t,x,y) 0 0 0.000 0 0.000 1/3 0.000 0.125 0.000 0.000 2/3 0.125 0.282 0.000 0.0430 1 0.407 0.563 0.0430 0.141 To numerically approximate y(1), where (1+t^2)y''+y'-y=0, y(0)=1,y'(0)=-1, using 4 steps of Euler's method, first convert to a system: y1' = y2, y1(0)=1; y2' = (y1-y2)/(1+t^2), y2(0)=-1. sage: RR = RealField(sci_not=0, prec=4, rnd='RNDU') sage: t, y1, y2=PolynomialRing(RR,3).gens() sage: f = y2; g = (y1-y2)/(1+t^2) sage: eulers_method_2x2(f,g, 0, 1, -1, 1/4, 1) t x h*f(t,x,y) y h*g(t,x,y) 0 1 -0.250 -1 0.500 1/4 0.750 -0.125 -0.500 0.282 1/2 0.625 -0.0546 -0.218 0.188 3/4 0.625 -0.00781 -0.0312 0.110 1 0.625 0.0196 0.0782 0.0704 To numerically approximate y(1), where y''+ty'+y=0, y(0)=1,y'(0)=0: sage: t,x,y=PolynomialRing(RR,3).gens() sage: f = y2; g = -y1-y2*t sage: eulers_method_2x2(f,g, 0, 1, 0, 1/4, 1) t x h*f(t,x,y) y h*g(t,x,y) 0 1 0.000 0 -0.250 1/4 1.00 -0.0625 -0.250 -0.234 1/2 0.938 -0.117 -0.468 -0.171 3/4 0.875 -0.156 -0.625 -0.101 1 0.750 -0.171 -0.687 -0.0156 AUTHOR: David Joyner """ print "%10s %20s %25s %20s %20s"%("t", "x","h*f(t,x,y)","y", "h*g(t,x,y)") n=int((RealField(max(6,RR.precision()))('1.0'))*(t1-t0)/h) t00 = t0; x00 = x0; y00 = y0 for i in range(n+ZZ(1)): print "%10r %20r %25r %20r %20r"%(t00,x00,h*f(t00,x00,y00),y00,h*g(t00,x00,y00)) x01 = x00 + h*f(t00,x00,y00) y00 = y00 + h*g(t00,x00,y00) x00 = x01 t00 = t00 + h