Sample Report On Numerical Solutions To Odes

The task is to solve ordinary differential equation using numerical methods.

The equation is:dydx=2xy+x3 (1)

Equation (1) is linear non-ordinary differential equation, let`s solve it.

First we need to solve the linear ordinary part, let`s rewrite it:
dydx=2xy

Using separation variables, we see:

dyy=2dxx
lny=2lnx+ lnC

Solve for y by exponentiating, and we get the general solution

y=x2C
As we found the general solution of linear ordinary part let`s find solution for non-oridnary part, to do this we should rewrite the general solution in next way
y=x2φx (2)
Where φx is an arbitrary function. Now we have to substitute (2) into (1), to find the unknown function φx.
2x φx+x2 φ'x=2xx2 φx+x3 (3)

Solving 3, we get that:

φx=x22+A

So the general solution of (1) is:

y=x2A+x42 (4)
2Ax+2x3=2xAx2+x42+x3

After simple operations we get that left part is equal to right part, so the general solution (4) is correct.

Using the initial condition y0.1=1 we see that A = 99.995, so the actual solution is:
y=12 (199.99x2+x4)

Let`s solve equation (1) using numerical methods.

The following code is written in Python:
import math
def euler(f, y0, a, b, h):
t, y = a, y0
result = []
while t < b + h:
result.append(y)
y += h * f(t, y)
t += h
return result
def midpoint(f, y0, a, b, h):
t, y = a, y0
result = []
while t < b + h:
result.append(y)
y += h * f(t + h / 2, y + h / 2 * f(t, y) )
t += h
return result
def rk4(f, y0, a, b, h):
t, y = a, y0
result = []
while t < b + h:
result.append(y)
k1 = h * f(t, y)
k2 = h * f(t + h / 2, y + k1 / 2)
k3 = h * f(t + h / 2, y + k2 / 2)
k4 = h * f(t + h, y + k3)
y += (k1 + 2 * k2 + 2 * k3 + k4) / 6
t += h
return result
def actual(f, y0, a, b, h):
t, y = a, y0
result = []
while t < b + h:
result.append(y)
t += h
y = (t * t * 199.99 + t ** 4) / 2
return result
def function(time, temp):
return (2 / time * temp + time ** 3)
h = 0.1
a = 0 + h
b = 2
y0 = 1
results = [ [], [], [], [] ]
results[0] = euler(function, y0, a, b, h)
results[1] = midpoint(function, y0, a, b, h)
results[2] = rk4(function, y0, a, b, h)
results[3] = actual(function, y0, a, b, h)
f = open('output.txt', 'w+')
t, y = a, y0
print "x\tEuler\t\tMidpoint\tRunge-Kutta-4\tActual"
for i in xrange(0, len(results[0])):
temp = "%2.1f\t%2.8f\t%2.8f\t%2.8f\t%2.8f" % (t, results[0][i], results[1][i], results[2][i], results[3][i])
print temp
f.write(temp + '\n')
t += h
print
print "Euler\t\tMidpoint\tRunge-Kutta-4"
for i in xrange(0, len(results[0])):
temp = "%2.8f\t%2.8f\t%2.8f" % (abs(results[0][i] - results[3][i]), abs(results[1][i] - results[3][i]), abs(results[2][i] - results[3][i]))
print temp
f.write(temp + '\n')
t += h
r = raw_input()

The results we can view at Figure 1, and the errors in Figure 2

Figure 1
Figure 2
As we can see the Euler method gives the worst performance, and the RK-4 method gives the best performance, because it`s error is the smallest.
All results were checked using wolfram alpha online tool.