Solving Ordinary Differential Equations (ODE)
- Jocelyne Zofia
- Jan 3, 2022
- 5 min read
Theoretical Physics - Computational Project
University of Fribourg, October 2012
AIM
Our program consists of 3 functions (corresponding to the 3 implemented methods), whose purpose is to find the solution of an ODE .
The 1st function (= odeeuler) uses an algorithm based on the Euler method in order to calculate the solution the Second function (= odemidp) is the implementation of the Midpoint method, whereas the third one, oderk4, is based on the Runge-Kutta 4th order method RK4.
Project report:
Code examples
from pylab import *
import numpy as np
import matplotlib.pyplot as plt
import math
def oderk4 (fxy, x0, y0, h, N):
x = x0;
y = y0;
k1= 0;
k2= 0;
k3= 0;
k4= 0;
E = 0.0;
Err = 0.0;
sol =[]
for i in range(1, N+1):
k1 = h * fxy(x, y)
k2 = h * fxy(x+0.5*h, y+0.5*k1)
k3 = h * fxy(x+0.5*h, y+0.5*k2)
k4 = h * fxy(x+h, y+k3)
y = y + (k1 + 2*k2 + 2*k3 +k4)/6
x += h
E += math.fabs(y-(math.exp(1-math.cos(x))))
Err = E/N
sol.append(y)
return sol, y
if __name__ == "__main__":
crit = 1
x0 = 0
y0 = 300
x1 = 10
N = 100
h = (x1-x0)/float(N)
yexp = 100
z0 = 0
while (crit > 0.01):
def fxy( x, y):
return sqrt(10)-(x*x)
x0 = 0
oderk4(fxy, x0, y0, h, N)
x0 = oderk4(fxy, x0, y0, h, N)[1]
def fxy( x, y):
return x
y0 = 0
oderk4(fxy, x0, y0, h, N)
Y = oderk4(fxy, x0, y0, h, N)[1]
crit = math.fabs(Y - yexp)
if Y < yexp:
z0 += 0.01
else:
z0 -= 0.005
y0 = 300+z0
print(Y, crit)
def fxy( x, y):
return x
y0 = 0
oderk4(fxy, x0, y0, h, N)
x = np.linspace(0, 10, N)
y = oderk4(fxy, x0, y0, h, N)[0]
plt.plot (x,y)
plt.xlabel('x')
plt.ylabel('y')
plt.savefig('DE26.png')
plt.show()
from pylab import *
import numpy as np
import matplotlib.pyplot as plt
import math
#The prodecure will be the same as in the previous progrms
'''
def odeeuler(fxy, x0, y0, h, N):
x = x0;
y = y0;
E = 0.0;
Err = 0.0;
sol =[]
for i in range(1, N+1):
y += h * fxy( x, y)
x += h
E += math.fabs(y-(math.exp(1-math.cos(x))))
Err = E/N
sol.append(y)
#return y
return sol, Err
'''
# We define the variables of the function (the local variables)
def oderk4 (fxy, x0, y0, h, N):
x = x0;
y = y0;
k1= 0;
k2= 0;
k3= 0;
k4= 0;
E = 0.0;
Err = 0.0;
#the function (object, to which I give the input)where the RK4 method is implemented
sol =[] #defines the 1-dimentional empty array
for i in range(1, N+1):
k1 = h * fxy(x, y)
k2 = h * fxy(x+0.5*h, y+0.5*k1)
k3 = h * fxy(x+0.5*h, y+0.5*k2)
k4 = h * fxy(x+h, y+k3)
y = y + (k1 + 2*k2 + 2*k3 +k4)/6
x += h
E += math.fabs(y-(math.exp(1-math.cos(x))))
Err = E/N
sol.append(y) #array in which I will have the stored solution 'y'.
#append means that it will add new result to the array
return sol, y #command giving the output of the difined function oderk4 (it will give me the array).
if __name__ == "__main__": #defines the main body
crit = 1 #deifnes the initial variables (global) (initializing of the value)
x0 = 0
y0 = 290 #we have a first shot of the initial condition z(0)
x1 = 10
N = 100
h = (x1-x0)/float(N)
yexp = 100
z0 = 0
#We define the presition of the calculation via the variable criterium
while (crit > 0.01):
def fxy( x, y):
return sqrt(10)-(x*x) #the second separated o.d.equation
x0 = 0
oderk4(fxy, x0, y0, h, N)
x0 = oderk4(fxy, x0, y0, h, N)[1]
def fxy( x, y): #define the first separated ODE dx/dx = z
return x # we defined z as x (according to our initial definition of the function oderk4, where we use what we previously defined)
y0 = 0
oderk4(fxy, x0, y0, h, N)
Y = oderk4(fxy, x0, y0, h, N)[1] #here we define the value Y as the actual solution of our ODE dx/dy
crit = math.fabs(Y - yexp) #we caluclate the new value of the criterium
#(criterium stands for (in our case) the distance from the expected value of the analytical solution
#A try of implementation of the bissection method
if Y < yexp: #numerical approach to the analytical expected value from the left
z0 += 0.01
else:
z0 -= 0.005 #same but from the right
y0 = 290+z0 #shooting of the actual z(0)
print(Y, crit) #visualisation of the 'criterium' after each step (to check if the program runs properly)
def fxy( x, y):
return x
y0 = 0
oderk4(fxy, x0, y0, h, N) #if the proper z(0) is found in the previous step, this fucntion calculates then the solution
#Now we are going to plot the solution
x = np.linspace(0, 10, N)
y = oderk4(fxy, x0, y0, h, N)[0]
plt.plot (x,y)
plt.xlabel('x')
plt.ylabel('y')
plt.savefig('DE26.png')
plt.show()
'''
def fxy( x, y):
return x
x0 = 0
y0 = 0
x1 = 10
N = 100
h = (x1-x0)/float(N)
crit = 1
yexp = 100
while (crit > 0.01):
oderk4(fxy, x0, y0, h, N)
x0 = y0
#y0 = 0
C = oderk4(fxy, x0, y0, h, N)[1]
crit = math.fabs(C - yexp)
print(C, crit)
#if oderk4(fxy, x0, y0, h, N)[1] < yexp:
def fxy(x, y):
return sqrt(100)-(x*x)
x0 = 0
#y0 = 0
x1 = 10
N = 100
h = (x1-x0)/float(N)
oderk4(fxy, x0, y0, h, N)[1]
if C < yexp:
y0 += yexp
#x0 = y0
else:
y0 -= y0/(2*N)
#x0 = y0
#print(C, crit)
'''
'''
def fxy(x, y):
return sqrt(C)-(x*x)
x0 = 0
y0 = 0
x1 = 10
N = 100
h = (x1-x0)/float(N)
yexp = 10
#crit = 1
K = 0
#while (crit > 0.01):
R = oderk4(fxy, x0, y0, h, N)[1]
SOL = oderk4(fxy, x0, y0, h, N)[0]
crit = math.fabs(oderk4(fxy, x0, y0, h, N)[1] - yexp)
K += K
if oderk4(fxy, x0, y0, h, N)[1] < yexp:
y0 += yexp
else:
y0 -= y0/(2*N)
print (R, crit)
'''
'''
yexp = 10
crit = 1
while (crit > 0.01):
def fxy( x, y):
return x
x0 = 0
y0 = 0
x1 = 10
N = 100
h = (x1-x0)/float(N)
C = oderk4(fxy, x0, y0, h, N)
def fxy(x, y):
return sqrt(C)-(x*x)
x1 = 10
N = 100
h = (x1-x0)/float(N)
oderk4(fxy, x0, y0, h, N)
crit = oderk4(fxy, x0, y0, h, N) - yexp
N += 1
if oderk4(fxy, x0, y0, h, N) < yexp:
y0 += yexp
else:
y0 -= y0/(2*N)
print (C)
#odeeuler(fxy, x0, y0, h, N)
#print (odeeuler(fxy, x0, y0, h, N))
#oderk4 (fxy, x0, y0, h, N)
#print(oderk4 (fxy, x0, y0, h, N))
#C = oderk4(fxy, x0, y0, h, N)
'''
'''
x = np.linspace(x0, x1, N)
y = SOL
plt.plot (x,y)
#plt.title ('Solution of y\' = y * sin(x) by Euler\'s method\n Division = '+str(h)+', Error =' +str(odeeuler(fxy, x0, y0, h, N)[1]))
plt.xlabel('x')
plt.ylabel('y')
plt.savefig('DE23.png')
plt.show()
'''
Comments