Week 8: Solving Differentital Equations

First lets take a simple relation ship, $$ \frac{dx}{dt} = -\alpha x $$

To solve this all we would have to do is realize that this is exponential decay and the solution is $$x(t) = e^{-\alpha t} $$

import numpy as np
import matplotlib.pyplot as plt
%config InlineBackend.figure_format='retina'

X=np.arange(0,10,0.01)

plt.plot(X,np.exp(-X))

plt.show()

def xfunc(x):
    return -x

We want to build a step in the derivative to approximate the exact solution

X=np.arange(0,2,0.01)

stepsize =0.4

x,y0=min(X),1

XX = [x]
YY = [y0]

while x< max(X):
    y_new = xfunc(YY[-1])*stepsize + YY[-1]
    YY.append(y_new)
    x+=stepsize
    XX.append(x)



plt.plot(X,np.exp(-X))
plt.plot(XX,YY,'.-','r')

plt.show()

X=np.arange(0,2,0.01)

for stepsize in np.linspace(0.05,0.5,5):
    
    x,y0=min(X),1
    XX = [x]
    YY = [y0]

    while x< max(X):
        y_new = xfunc(YY[-1])*stepsize + YY[-1]
        YY.append(y_new)
        x+=stepsize
        XX.append(x)
    plt.plot(XX,YY,'.-',label = str(round(stepsize,2)), alpha = 0.3)
    


plt.plot(X,np.exp(-X), label='exact solution')

plt.legend()

plt.show()

Lets Generalize this approach

$$\frac{dx}{dt}=f(x,t)$$

We want $f(x,t)$ to be dependent on position and time.

We want to plot a solution of this first order differential equation

We have the following inputs

we have a function $f(x,t)$

inital position at time zero

a stepsize

End time

def xfunc(x,t):
    return -x

def ode_1d_solve(xfunc,X0,dt,T):
    Num = int(T//dt)
    x = np.ones(Num +1)*X0
    t = np.linspace(0,Num*dt,len(x))
    for count in range(Num):
        x[count+1]=x[count]+xfunc(x[count],t[count])*dt
        
    return t,x

t,x = ode_1d_solve(xfunc,1,0.2,2)

plt.plot(t,x,'.-')
plt.show()

%%time
X=np.arange(0,2,0.01)

t,x = ode_1d_solve(xfunc,1,0.02,2)

plt.plot(t,x,'.','r')


plt.plot(X,np.exp(-X), 'b',label='exact solution')

plt.legend()

plt.show()

CPU times: user 200 ms, sys: 15.4 ms, total: 216 ms
Wall time: 242 ms

len(t)

100000

%%time

def xfunc(x,t):
    return (np.exp(t/3)-x)/2



t,x = ode_1d_solve(xfunc,1,0.02,20)

plt.plot(t,x,'r','.')

plt.show()

CPU times: user 169 ms, sys: 13.1 ms, total: 182 ms
Wall time: 198 ms

Lets look at the logistic growth, in the continuous model

$$ \frac{dx}{dt}=rx(K-x)$$

def logistic_func(x,t):
    return 0.15*x*(10-x)

t,x = ode_1d_solve(logistic_func,1,0.01,5)

plt.plot(t,x,'r')

plt.show()

plt.figure(figsize=(13,5))

for s in np.linspace(1,21,11):
    t,x = ode_1d_solve(logistic_func,s,0.001,6)
    
    plt.plot(t,x, label='starting value = '+str(round(s,1)))
    
plt.legend(loc='center left', bbox_to_anchor=(1,0.5))
plt.show()