2

I tried to write an algorithm for solving the non-linear ODE 

dr/dt = r(t)+r²(t)

which has (one possible) solution 

r(t) = r²(t)/2+r³(t)/3

For that I implemented both the euler method and the RK4 method in python. For error checking I used the example on rosettacode: 

dT/dt = -k(T(t)-T0)

with the solution 

T0 + (Ts − T0)*exp(−kt)

Thus, my code now looks like

import numpy as np
from matplotlib import pyplot as plt

def test_func_1(t, x):
    return x*x

def test_func_1_sol(t, x):
    return x*x*x/3.0

def test_func_2_sol(TR, T0, k, t):
    return TR + (T0-TR)*np.exp(-0.07*t)

def rk4(func, dh, x0, t0):
    k1 = dh*func(t0, x0)
    k2 = dh*func(t0+dh*0.5, x0+0.5*k1)
    k3 = dh*func(t0+dh*0.5, x0+0.5*k2)
    k4 = dh*func(t0+dh, x0+k3)
    return x0+k1/6.0+k2/3.0+k3/3.0+k4/6.0

def euler(func, x0, t0, dh):
    return x0 + dh*func(t0, x0)

def rho_test(t0, rho0):
    return rho0 + rho0*rho0

def rho_sol(t0, rho0):
    return rho0*rho0*rho0/3.0+rho0*rho0/2.0

def euler2(f,y0,a,b,h):
    t,y = a,y0
    while t <= b:
        #print "%6.3f %6.5f" % (t,y)
        t += h
        y += h * f(t,y)

def newtoncooling(time, temp):
    return -0.07 * (temp - 20)

x0 = 100
x_vec_rk = []
x_vec_euler = []
x_vec_rk.append(x0)

h = 1e-3
for i in range(100000):
    x0 = rk4(newtoncooling, h, x0, i*h)
    x_vec_rk.append(x0)

x0 = 100
x_vec_euler.append(x0)
x_vec_sol = []
x_vec_sol.append(x0)

for i in range(100000):
    x0 = euler(newtoncooling, x0, 0, h)
    #print(i, x0)
    x_vec_euler.append(x0)
    x_vec_sol.append(test_func_2_sol(20, 100, 0, i*h))

euler2(newtoncooling, 0, 0, 1, 1e-4)

x_vec = np.linspace(0, 1, len(x_vec_euler))

plt.plot(x_vec, x_vec_euler, x_vec, x_vec_sol, x_vec, x_vec_rk)
plt.show()

#rho-function
x0 = 1
x_vec_rk = []
x_vec_euler = []
x_vec_rk.append(x0)

h = 1e-3
num_steps = 650
for i in range(num_steps):
    x0 = rk4(rho_test, h, x0, i*h)
    print "%6.3f %6.5f" % (i*h, x0)
    x_vec_rk.append(x0)

x0 = 1
x_vec_euler.append(x0)
x_vec_sol = []
x_vec_sol.append(x0)

for i in range(num_steps):
    x0 = euler(rho_test, x0, 0, h)
    print "%6.3f %6.5f" % (i*h, x0)
    x_vec_euler.append(x0)
    x_vec_sol.append(rho_sol(i*h, i*h+x0))

x_vec = np.linspace(0, num_steps*h, len(x_vec_euler))
plt.plot(x_vec, x_vec_euler, x_vec, x_vec_sol, x_vec, x_vec_rk)
plt.show()

It works fine for the example from rosettacode, but it is unstable and explodes (for t > 0.65, both for RK4 and euler) for my formula. Is therefore my implementation incorrect, or is there another error I do not see?

arc_lupus
  • 3,942
  • 5
  • 45
  • 81

1 Answers1

2

Searching for the exact solution of your equation:

dr/dt = r(t)+r²(t)

I've found:

r(t) = exp(C+t)/(1-exp(C+t))

Where C is an arbitrary constant which depends on initial conditions. It can be seen that for t -> -C the r(t) -> infinity.

I do not know what initial condition you use, but may be you meet with this singularity when calculate numerical solution.

UPDATE: Since your initial condition is r(0)=1 the constant C is C = ln(1/2) ~ -0.693. It can explain why your numerical solution crashes at some t>0.65

UPDATE: To verify your code you can simply compare your numerical solution calculated for, say 0<=t<0.6 with exact solution.

Roman Fursenko
  • 688
  • 4
  • 9
  • My initial condition is x0 = 1, as seen in the code. – arc_lupus Feb 16 '17 at 11:34
  • You can work in the initial value to get `r(t) = ( r(0)*exp(t) )/( 1-r(0)*(exp(t)-1) )`. [WA](http://www.wolframalpha.com/input/?i=r%27%28t%29+%3D+r%28t%29^2%2Br%28t%29,+r%280%29%3Da) – Lutz Lehmann Feb 16 '17 at 11:38
  • @arc_lupus, Also, to verify your code you can simply compare your numerical solution calculated for, say `0<=t<0.6` with exact solution. – Roman Fursenko Feb 16 '17 at 11:55