0
$\begingroup$

This is my first post on here, so please excuse mistakes if any.

I am trying to plot out the difference between two ODE solvers based on Taylor series:

  • 1st order acccurate: $x(t_0 + h) = x(t_0) + h x'(t_0)$ and
  • 2nd order accurate: $x(t_0 + h) = x(t_0) + h x'(t_0) + \frac{h^2}{2} x''(t_0)$

ODE: x'= -x + e^-t x0 = 0; t0 = 0 x'' = x - 2e^-t

For some reason real solution and both 1st and 2nd order approximation are the same. Any ideas?

from math import exp
from numpy import arange
from numpy import zeros
from matplotlib import pyplot as plt

def f(t,x):
    dx = -x + exp(-t)
    return dx

def df(t,x):
    df = x - 2*exp(-t)
    return df

def real(h):
    t = 0
    xk_list = []
    tk_list = []
    def func(t):
        x = t*exp(-t)
        return x

    while(True):
        xk_list.append(func(t))
        tk_list.append(t)
        t += h
        if (t > 3.0):
            return tk_list, xk_list

def forward_euler_1(x0, t0, h):
    xk = x0
    tk = t0
    tk_list = []
    xk_list = []
    counter = 0
    while(True):
        xk1 = xk + h*f(tk,xk)
        xk_list.append(xk)
        tk_list.append(tk)
        if(tk > 3.0):
            print(counter,tk)
            return tk_list, xk_list
        xk = xk1
        tk = tk + h
        counter += 1

def forward_euler_2(x0,t0,h):
    xk = x0
    tk = t0
    tk_list = []
    xk_list = []
    counter = 0 

    while(True):
        xk1 = xk + h*f(tk, xk) + (h**2)/2*df(tk, xk)
        xk_list.append(xk1)
        tk_list.append(tk)    
        if(tk > 3.0):
            print(counter,tk)
            return tk_list, xk_list        
        xk = xk1
        tk = tk + h
        counter += 1

x, y = forward_euler_1(0.0, 0.0,0.001)
x2, y2 = forward_euler_2(0.0, 0.0, 0.001)
x3, y3 = real(0.001)

plt.plot(x,y)
plt.plot(x2,y2,'--')
plt.plot(x3, y3)
plt.show()       
$\endgroup$
  • $\begingroup$ I'm an idiot, there is nothing wrong with the code, I am just using an h resoution that is too high, if one put's the time discretization at 0.1 you see what I expected to see. $\endgroup$ – terraregina Jan 19 '19 at 13:03
  • 1
    $\begingroup$ Hi welcome to scicomp.stackexchange. If your comment answers (and resolve) your question you can post an answer and accept it. $\endgroup$ – Mauro Vanzetto Jan 19 '19 at 14:37
  • $\begingroup$ But delete the "I'm an idiot" part first ;-) Welcome to scicomp SE! $\endgroup$ – GertVdE Jan 22 '19 at 20:18
  • $\begingroup$ done and thank you very much, it has been a most gracious welcome ;D $\endgroup$ – terraregina Jan 23 '19 at 14:36
0
$\begingroup$

There is nothing wrong with the code, I am just using an h resoution that is too high, if one put's the time discretization at 0.1 you see what I expected to see.

| cite | improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.