Solving ODE with Taylor method gives very wrong approximation

Question

Can someone check this python code. I need to use Taylor's method of order 2 to approximate the solution to $$ y'= \frac1{x^2}-\frac{y}{x}-y^2,~~ 1\le x\le 2,~~ y(1)=-1 ~\text{ and }~ h=0.05. $$ it gives me very big approximate number and a wrong sign. The exact is $y(x)= -1/x$, when $x=1.1$, $y=-9.090909091$

# Python Code to find the approximation of an ordinary
# differential equation using Taylor method.

# Given
# dy / dx =(1/x^2)-(y/x)-(y^2), y(1)=-1, h=0.05
def func(x, y):
    return (1/(x**2))-(y/x)-(y**2)


# Function for euler formula
def euler(x0, y, h, x):
    temp = -0



    # Iterating till the point at which we
    # need approximation
    while x0 < x:
        temp = y
        y = (1/(x**2))-(y/x)-(y**2)
        x0 = x0 + h

x0 = 1
y0 = -1
h = 0.05

# Value of x at which we need approximation
x = 1.1

euler(x0, y0, h, x)


temp=-0
def second_order(x0,y,h,x):
    while x0 < x:
        temp = y
        y = (3/(x**3))+(3*(y**2)/x)+2*(y**3)
        x0 = x0 + h
    print("Approximate solution at x = ", x, " is ", "%.6f" % y)

second_order(x0,y0,h,x)

Answer 1

Removing code that isn't used or doesn't do anything, you have

x0 = 1
y0 = -1
h = 0.05

# Value of x at which we need approximation
x = 1.1

def second_order(x0,y,h,x):
    while x0 < x:
        y = (3/(x**3))+(3*(y**2)/x)+2*(y**3)
        x0 = x0 + h
    print("Approximate solution at x = ", x, " is ", "%.6f" % y)

second_order(x0,y0,h,x)

This code does not compute $y(x)$ as x0 walks from $1$ to $1.1$. (I wish there were an intelligible way to write that sentence. Here, let me rewrite your code to do what it appears it wants to do but so that the variable names actually correspond to the semantics of their values.)

xStart = 1
yStart = -1
xEnd = 1.1
h = 0.05

x = xStart
y = yStart
while x < xEnd:  # This only allows the last 
    # pass through the loop because the 
    # internal represenation of 0.05 is very 
    # slightly less than 1/20.  Probably 
    # better to use x <= xEnd if you want x to 
    # *reach* xEnd for various values of 
    # xStart, xEnd, and h.
    y = (3/(x**3))+(3*(y**2)/x)+2*(y**3)
    x = x + h
    #Uncomment the next line to get your own table of intermediates.
    #print(x,y)
print("Approximate solution at x = ", x, " is ", "%.6f" % y)

Now we have $(x,y)$ tracing out a curve where xStart${} < x < {}$xEnd. Which curve? Well, let $x_0 = {}$xStart${} = 1$ and $y_0 = {}$yStart${} = -1$. For $1 \leq i \leq 2$, let $x_i = x_0 + i$h${} = x_0 + i/20$ and $y_i = \frac{3}{x_{i-1}^3} + \frac{3 y_{i-1}^2}{x_{i-1}} + 2 y_{i-1}^3$.

Let's make a table of values.

\begin{align*} &i & &x_i & &y_i \\ \hline &0 & &1 & &-1 \\ &1 & &1.05 & &4 \\ &2 & &1.10 & &\frac{544\,256}{3087} = 176.306{\dots} \end{align*} So this Python code will give you a positive number as result.

Searching the Internet for uses of "taylor method update numerical differential equation" and similar phrases turns up no hits, so it is unclear what method you are attempting to implement. (There are several reference to Taylor series, but nothing you write looks like you are using the Taylor series of $-1/x$ centered at $1$, which is $$ T_{-1/x}(x) = -1 + (x-1) - (x-1)^2 + (x-1)^3 - \cdots \text{,} $$ so those don't seem relevant.) Since I can't guess what variant of an Euler integrator you are intending to use (and there are many, many, MANY, MANY such variants), I cannot make recommendations to steer your code back to your intended calculation.

Answer 2

Let's first check the exact solution that is given as reference solution, sometimes it happens that is wrong, from another task or such. Another concern is that sometimes the IVP solution is unstable in that close-by solutions at the initial point develop qualitatively different, moving away very fast. After that I'll address the numerical issues.

---

Your equation $$y'= \frac1{x^2}-\frac{y}{x}-y^2$$ is a Riccati equation. Solving it via $y=\frac{u'}{u}$ results in the equation $$ x^2u''+xu'-u=0,~~ u(1)=1,~u'(1)=-1 $$ which is now an Euler-Cauchy equation with basis solutions of type $y=x^m$. Inserting gives $m=\pm 1$, so that the full solution is $$ u(x)=Ax+Bx^{-1}\implies y(x)=\frac{A-Bx^{-2}}{Ax+Bx^{-1}}=\frac{Ax^2-B}{x(Ax^2+B)} $$ which for $A=0$, $B=1$ results in the IVP solution. This solution formula appears to be stable under perturbation of the initial conditions, so that there should be no numerical surprises from this side.

---

The Taylor methods use the evaluation of the Taylor polynomial for $y(x+h)$ at $x$ to step forward. If the linear Taylor polynomial is used, the explicit Euler method results. Here it is demanded to use the quadratic Taylor polynomial. $$y(x+h)\approx y(x)+hy'(x)+\frac12h^2y''(x).$$ To apply this one needs to compute the second derivative of $y$. This can be done by differentiating the given ODE (using the chain rule) as per $$ y'(x)=f(x,y(x))\implies y''(x)=\partial_xf(x,y(x))+\partial_yf(x,y(x))f(x,y(x)). $$ This method is no longer Runge-Kutta, as more than the values of the ODE function are used.

---

The easiest way to implement this in a human-readable fashion is to compute the partial derivatives separately instead of inserting into one huge expression.

def f(x,y): return x**-2 - y/x - y**2 
def f_x(x,y): return -2*x**-3 + y/x**2 
def f_y(x,y): return - 1/x - 2*y
def Taylor2step(x,y,h):
    Dy = f(x,y)
    D2y = f_x(x,y)+f_y(x,y)Dy
    return y+hDy+0.5h2D2y

and then iterate over the time span

while x < xf+0.1*h:
    x,y = x+h, Taylor2step(x,y,h) 

In the first iterations this gives

x=  1.0000:  y= -1.000000000000
x=  1.0500:  y= -0.952500000000
x=  1.1000:  y= -0.909313789767

This looks correct for an error term $O(th^2)$.