# -*- coding: utf-8 -*-
"""
Created on Fri May 08 10:47:01 2015
@author: Kevi
"""
from __future__ import division
from scipy import *
import matplotlib.pyplot as plt
import numpy as np
from scipy import special
from scipy import optimize
from math import e
from scipy import integrate
from scipy.integrate import odeint
""" Hi and welcome to this IVP solver program. To be able to run this program
you have to have an IVP problem (or system of them) in the following form
y´=f(t,y) with y0=f(t0). This code is a 2-step explicit method. This code uses
splines of degree two to illustarte the desire result.
"""
""" The first thing you have to do to run this program is to write your
function. To solve your problem you need to write your IVP below in the
constructed code. Your function should be an array. The dimension of the array
will depend on how your function looks like. If you just have one equation to
solve write it in the follwing way "array([-1*y[0]])". The y[0] tells the program that to
call that specific y. If you have a system of equationfor example a 2x2 matrix
of the form
|-y y |
|y -3y|
then you write "array([-y[0] +1*y[1],y[0]-3*y[1]])"
Now you are half way to calculate your function"""
def function(y,t):
""" Give your IVP problem here."""
return array([5*y[0]])
class IVP_Solver(object):
""" To run the program and get the answear for your problem you have to
call the
program ´a=IVP_solver(function,y_condition,theta,TOL,perc)´ then a(t_start.t_end)
The program returns a plot of your answear in your interval and how many points
it took to get there.
"""
def __init__(self,function,y_condition,theta,TOL,perc):
self.y_condition=y_condition
self.function=function
self.theta=theta
self.TOL=TOL
self.perc=perc
def __call__(self,t_start,t_end):
return self.Explicit_2_step_method_with_splines(self.function,self.y_condition,t_start,t_end,self.theta,self.TOL,self.perc)
def error_const(self,theta,eta):
"""Calculates the error constant.
"""
self.theta=theta
self.eta=eta
self.K=0
self.k=len(theta) +1
if self.k==2:
self.cx= cos(theta)
self.sx= sin(theta)
self.cy = cos(eta)
self.sy = sin(eta)
self.Kx=(2*self.cx-5*self.sx)/(self.cx-2*self.sx)/6
self.Ky=(2*self.cy-5*self.sy)/(self.cy-2*self.sy)/6
self.K= self.Kx/(self.Ky-self.Kx)
return self.K
def pol(self,theta,h,y,yp):
"""Returns a points which helps with constructing a spline between two
given points
"""
self.h=h
self.y=y
self.yp=yp
self.kk= len(theta) +1
if self.kk ==2:
self.c=cos(theta)
self.s=sin(theta)
self.p = (self.c*(y[-2]-y[-1]+h[-2]*yp[-1])+self.s*h[-2]*(yp[-2]-yp[-1]))/(h[-2]**2*(self.c-2*self.s))
return self.p
def angle2(self,theta):
"""Returns eta, eta depends what your theta is. If theta=pi/2 the method
will be Adam Bashforth."""
self.theta=theta
self.k= len(theta)+1
self.eta=theta.copy()
for j in range(self.k-1):
if theta[j] <= 5*pi/8:
self.eta[j]=13*pi/16
else:
self.eta[j]=7*pi/16
return self.eta
def start_(self,theta,t,h,y,yp):
"""This function helps you to start with the 2-step explicit method.
Since in your problem you are just given you one point(we need two points to
start the process) this program calculates the second point with help of the
inbuild function "odeint" (Runge-Kutta). It also gives you the
second point in time you need to have."
"""
self.theta=theta
self.t=t
self.h=h
self.y=y
self.yp=yp
self.k = len(theta)+1
for j in range(self.k-1):
t.append(t[-1]+h[-1])
self.Y=odeint(self.function,y[0],t)
y.append(self.Y[-1])
self.f=self.function(y[-1],t[-1])
yp.append(self.f)
return t,y,yp
def controller(self,parvec,theta,k, TOL, est, r, perc,w):
"""The function controller choses how big the next stepsize should be.
If the given function grows fast or slow in an interval the controller will
take an smaller h, if the function does not change much in an interval the
controller will take a bigger h. This function estimates the error at easch
step, and tries to keep it equal to the tolerance"""
self.ord_ = len(theta) + 1
if self.k== self.ord_ +1:
self.cerrold=1
else:
self.cerrold= (TOL/est[-2])**(1/(self.ord_+1))
if est[-1]==0:
est[-1]=1e-18
self.cerr= (TOL/(est[-1]))**(1/(self.ord_+1))
self.convec =array([self.cerr,self.cerrold,r])
r= prod(array([self.convec[0]**self.parvec[0],self.convec[1]**self.parvec[1],self.convec[2]**self.parvec[2]]))
if r<perc:
r=perc
w=w+1
elif r > 2-perc:
r=2-perc
w=w+1
return r,w
def Explicit_2_step_method_with_splines(self,function,y_condition,t_start,t_end,theta,TOL,perc):
"""This program uses other programs below. The program returns a plot of
the numerical answear for your ODE problem and returns how many points it took to get there.
This process is an explicit method, the process stopiterating until the given
time interval is reached. Changing theta will give different explicit methods,
some are zero-stable some are not. "theta=array([pi/2]" is Adams-Bashforth method.
"""
self.function=function
self.y_condition=y_condition
self.t_start=t_start
self.t_end=t_end
self.parvec=array([1,1,0])/6
self.theta=array([theta])
self.ord_=len(self.theta)+1
self.TOL = TOL
self.perc=perc
self.w=0
self.eta=self.angle2(self.theta)
self.K=self.error_const(self.theta,self.eta)
self.h=[self.TOL**(2/(self.ord_+1))]
self.t=[t_start]
self.y=[y_condition]
self.yp=[self.function(y_condition,self.t[0])]
self.t,self.y,self.yp=self.start_(self.theta,self.t,self.h,self.y,self.yp)
self.est=[0]
self.r=1
self.k=self.ord_
while self.t[self.k-1]<t_end:
self.k=self.k+1
self.h.append(self.r*self.h[self.k-3])
self.t.append(self.t[self.k-2]+self.h[self.k-2])
if self.t[-1] >t_end:
self.t[-1]=t_end
self.px = self.pol(self.theta,self.h,self.y,self.yp)
self.py = self.pol(self.eta,self.h,self.y,self.yp)
self.y.append(self.px*self.h[self.k-2]**2 +self.yp[self.k-2]*self.h[self.k-2]+self.y[self.k-2])
self.yp.append(self.function(self.y[self.k-1],self.t[self.k-1]))
self.y_1= self.py*self.h[self.k-2]**2 +self.yp[self.k-2]*self.h[self.k-2]+self.y[self.k-2]
self.est.append(norm(self.K*(self.y_1-self.y[self.k-1])/(abs(self.y[self.k-1])+1e-3)))
self.r,self.w=self.controller(self.parvec,self.theta,self.k, self.TOL, self.est, self.r, self.perc,self.w)
if len(self.t) >=100000:
return "To big interval or wrong function!"
self.y_finale=zeros(shape=(len(self.t),len(self.y_condition)))
for j in range(len(self.t)):
for k in range(len(self.y_condition)):
self.y_finale[j,k]=self.y[j][k]
for i in range(len(self.y_condition)):
plt.plot(array(self.t),self.y_finale[:,i])
plt.ylabel('Y-VALUES')
plt.xlabel('TIME')
plt.title('Result of your IVP')
plt.show()
return len(self.t)