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 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 at 14:37
  • $\begingroup$ But delete the "I'm an idiot" part first ;-) Welcome to scicomp SE! $\endgroup$ – GertVdE Jan 22 at 20:18
  • $\begingroup$ done and thank you very much, it has been a most gracious welcome ;D $\endgroup$ – terraregina Jan 23 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.

$\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.