CODEWITHSUNDEEP
python-3.xnumerical-methodsode
user7918565
user7918565

Reputation: 11

Solving 2nd order ODE with python

I was writing some code to solve 2nd order differential equation, but it gives a completely wrong result. I think that the problem is in expressing Euler method in the right way. I also tried another 2nd order ODE, but I also failed at approximating y(x). Could you point where the mistake could be? Please see the graph and the code:

Solving ODE:

y"(x)=-1001y'(x)-1000y(x)

y(0)=1, y'(0)=0

Analytical solution:

y(x)=(1000exp(-x)-exp(-1000x))/999

Rewriting as a system of 2 ODEs;

y'(x)=v(x)

v'(x)=-1001v(x)-1000y(x)=f(x,y,v)

y(0)=1,v(0)=0

The code in Jupyter Notebook (Python):

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

a=-0.0020
b=0.0
h=0.0005 #step size
x=np.arange(a,b,h)
n=int((b-a)/h)

def y_exact(x):
    return (1000*np.exp(-x)-np.exp(-1000*x))/999
def f(y,v):
    return -1001*v-1000*y

y_eu=np.zeros(n,float)
y_ei=np.zeros(n,float)

y_eu[0]=1.0
y_ei[0]=1.0

v_eu[0]=0.0
v_ei[0]=0.0

for j in range(0,n-1):
    #Euler's method (explicit)
    y_eu[j+1]=y_eu[j]+v_eu[j]*h
    v_eu[j+1]=v_eu[j]+f(y_eu[j],v_eu[j])*h

    #Implicit Euler's method
    v_ei[j+1]=(v_ei[j]-1000*y_ei[j]*h)/(1+1001*h+1000*h**2)   
    y_ei[j+1]=y_ei[j]+v_ei[j+1]*h
 plt.plot(x,y_exact(x),label="Exact solution")
 plt.plot(x,y_eu,label="Euler's method")
 plt.plot(x,y_ei,label="Implicit Euler's method")
 plt.legend()
 plt.show()

Thanks!

Upvotes: 0

Views: 6083

Answers (1)

user7918565
user7918565

Reputation: 11

Thanks to commenter Maarten Fabre's answer. The code did not work properly because initial conditions, that is y0 and v0 for Euler's method were taken wrongly:

def y_exact(x):
    return (1000*np.exp(-x)-np.exp(-1000*x))/999
def dydx_exact(x):
    return -1000*(np.exp(-x)-np.exp(-1000*x))/999

y_eu[0]=y_exact(a)
v_eu[0]=dydx_exact(a)

Upvotes: 1

Related Questions