1

I have the necessity to check if two splines are close enough with a certain tolerance to each other. the first spline:

spl_psi=scipy.interpolate.UnivariateSpline(t,psi,k=4,s=0)

is the "data input" one used to check the second, which is:

def fr1(x,y,z):
   return z

def fr2(t,r1,r2):
   return ((K_1*spl_delta(t)+(K_1*T_3)*spl_delta1(t)-(T_1*T_2)*r2-r1)/(T_1+T_2))

for i in range(0,len(T2)-1,1):  
    k1=h*fr1(T2[i],r1[i],r2[i])
    l1=h*fr2(T2[i],r1[i],r2[i])
    k2=h*fr1(T2[i]+0.5*h,r1[i]+0.5*k1,r2[i]+0.5*l1)
    l2=h*fr2(T2[i]+0.5*h,r1[i]+0.5*k1,r2[i]+0.5*l1)
    k3=h*fr1(T2[i]+0.5*h,r1[i]+0.5*k2,r2[i]+0.5*l2)
    l3=h*fr2(T2[i]+0.5*h,r1[i]+0.5*k2,r2[i]+0.5*l2)
    k4=h*fr1(T2[i]+h,r1[i]+k3,r2[i]+l3)
    l4=h*fr2(T2[i]+h,r1[i]+k3,r2[i]+l3)
    r1[i+1]=r1[i]+(k1+2*k2+2*k3+k4)/6
    r2[i+1]=r2[i]+(l1+2*l2+2*l3+l4)/6

spl_r1_n2=scipy.interpolate.UnivariateSpline(T2,r1,k=5,s=0) 
spl_psi_n2=scipy.interpolate.UnivariateSpline.antiderivative(spl_r1_n2)

spl_psi_n2 is my second spline that must be close to spl_psi, and T_1,T_2,T_3,K_1 are my parameters. This method is called Runge-Kutta and it's used to solve differential equation such as mine. Now, I want to find the parameters that make the best spline and then save them in an array called "res". This is my code: 

def fr1(x,y,z):
return z

def fr2(x,y,z):
return ((K*spl_delta(x)+(K*(t_3/4))*spl_delta1(x)-((t_1/4)+(t_2/4))*z-y)/(t_1*t_2/(16) ))   

r1_0=spl_psi1(0)
r2_0=spl_psi2(0)
h=T2[2]-T2[1]
r1=np.array( [r1_0]*len(T2))
r2=np.array( [r2_0]*len(T2))
toll=0.5
res=[]

for t_1 in range (2,200,1):
  for t_2 in range (2,50,1):
    for t_3 in range (2,50,1):
        for i in range(0,len(T2)-1,1):
            k1=h*fr1(T2[i],r1[i],r2[i])
            l1=h*fr2(T2[i],r1[i],r2[i])
            k2=h*fr1(T2[i]+0.5*h,r1[i]+0.5*k1,r2[i]+0.5*l1)
            l2=h*fr2(T2[i]+0.5*h,r1[i]+0.5*k1,r2[i]+0.5*l1)
            k3=h*fr1(T2[i]+0.5*h,r1[i]+0.5*k2,r2[i]+0.5*l2)
            l3=h*fr2(T2[i]+0.5*h,r1[i]+0.5*k2,r2[i]+0.5*l2)
            k4=h*fr1(T2[i]+h,r1[i]+k3,r2[i]+l3)
            l4=h*fr2(T2[i]+h,r1[i]+k3,r2[i]+l3)
            r1[i+1]=r1[i]+(k1+2*k2+2*k3+k4)/6
            r2[i+1]=r2[i]+(l1+2*l2+2*l3+l4)/6
        spl_r1_n2=scipy.interpolate.UnivariateSpline(T2,r1,k=5,s=0)
        spl_psi_n2=scipy.interpolate.UnivariateSpline.antiderivative(spl_r1_n2)
        print ([t_1,t_2,t_3])
        for i in range (0,len(T2)-1,1):
            if fabs(spl_psi_n2(T2[i])-spl_psi(T2[i]))<toll:
                print 'ok'
                continue
            elif fabs(spl_psi_n2(T2[i])-spl_psi(T2[i]))>toll:
                print 'nope',i
                break
            elif i==len(T2)-2:
                print 'it's done',t_1,t_2,t_3
                res.append([t_1,t_2,t_3])
                break

Now, my issue is that when I run this it gives me 'It's done' for every combination of t_1,t_2,t_3. I am quite sure that the problem is in the chunk: 

for i in range (0,len(T2)-1,1):
    if fabs(spl_psi_n2(T2[i])-spl_psi(T2[i]))<toll:
        print 'ok'
        continue
    elif fabs(spl_psi_n2(T2[i])-spl_psi(T2[i]))>toll:
        print 'nope',i
        break
    elif i==len(T2)-2:
        print 'it's done',t_1,t_2,t_3
        res.append([t_1,t_2,t_3])
        break

but can't spin my way around it. Any help?

Enrico Pascatti
  • 31
  • 4

0 Answers0