# -*- 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 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.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)