value_with_derivatives

def value_with_derivatives (   self,   x   ) get the value of the function, and its first & second derivative Parameters: x  the point at which to evaluate the function Returns: (f(x), f'(x), f''(x)) Reimplemented from C2Ratio. Definition at line 1569 of file C2Functions.py. 01569 : 01570 "combine left and right functions into ratio, being very careful about zeros of the denominator" 01571 cache=self.cache 01572 if cache is None or x < cache[0] or x > cache[2]: #can't get it out of cache, must compute something 01573 y0, yp0, ypp0=self.left.value_with_derivatives(x) 01574 y1, yp1, ypp1=self.right.value_with_derivatives(x) 01575 01576 if ( abs(y1) < self.eps0*abs(yp1)*(abs(x)+self.dblmin) ): # paranoid check for relative proximity to a simple denominator zero 01577 dx0=self.dxsolve(x, y0, yp0, ypp0) 01578 dx1=self.dxsolve(x, y1, yp1, ypp1) 01579 01580 if abs(dx1-dx0) > self.eps1*abs(x): 01581 raise C2Exception("y/0 found at x=%g in LHopitalRatio" % x) 01582 01583 dx=abs(x)*self.epsx1+self.epsx2 01584 01585 #compute the walls of the region to be fixed with l'Hopital 01586 x1=x+dx1 01587 x0=x1-dx 01588 x2=x1+dx 01589 01590 y0, yp0, ypp0=C2Ratio.value_with_derivatives(self, x0) #get left wall values from conventional computation 01591 y2, yp2, ypp2=C2Ratio.value_with_derivatives(self, x2) #get right wall values 01592 01593 yy0, yyp0, yypp0=self.left.value_with_derivatives(x1) 01594 yy1, yyp1, yypp1=self.right.value_with_derivatives(x1) 01595 #now use L'Hopital's rule to find function at the center 01596 y1=yyp0/yyp1 01597 01598 #scale derivs to put function on [-1,1] since mma solution is done this way 01599 yp0*=dx 01600 yp2*=dx 01601 ypp0*=dx*dx 01602 ypp2*=dx*dx 01603 01604 #y[x_] = y1 + x (a + b x) + (x-1) x (x+1) (c + d x + e x^2 + f x^3) 01605 coefs=( y1, 01606 -(y0 - y2)/2., 01607 (y0 - 2*y1 + y2)/2., 01608 (7*y0 - 7*y2 + 7*yp0 + 7*yp2 + ypp0 - ypp2)/16., 01609 (-16*y0 + 32*y1 - 16*y2 - 9*yp0 + 9*yp2 - ypp0 - ypp2)/16., 01610 (-3*y0 + 3*y2 - 3*yp0 - 3*yp2 - ypp0 + ypp2)/16., 01611 (8*y0 - 16*y1 + 8*y2 + 5*yp0 - 5*yp2 + ypp0 + ypp2)/16. 01612 ) 01613 #y'[x] = a + 2 b x + (3x^2 - 1) (c + d x + e x^2 + f x^3) + (x-1) x (x+1) (d + 2 e x + 3 f x^2 ) 01614 #y''[x] = 2b + (x-1) x (x+1) (2 e + 6 f x) + 2 (3 x^2 -1) (d + 2 e x + 3 f x^2 ) + 6 x (c + d x + e x^2 + f x^3) 01615 01616 self.cache=x0,x1,x2, coefs 01617 01618 else: #not close to a zero of the denominator... compute conventional answer for ratio 01619 return y0/y1, (yp0*y1-y0*yp1)/(y1*y1), (y1*y1*ypp0+y0*(2*yp1*yp1-y1*ypp1)-2*y1*yp0*yp1)/(y1*y1*y1) 01620 01621 #if we get here, the poly coefficients are ready to go. 01622 x0, x1, x2, (y1, a,b,c,d,e,f)=self.cache 01623 01624 dx0=x1-x0 01625 dx=(x-x1)/dx0 01626 01627 q1=c + dx*(d + dx*(e + dx*f)) 01628 q2 =d + dx*(2*e + dx*3*f) 01629 q3=2*e+6*f*dx 01630 01631 xp1=(dx-1)*(dx+1)*dx 01632 xp2=(3*dx*dx-1) 01633 01634 y= y1 + dx*(a+b*dx) + xp1*q1 01635 yp=a + 2*b*dx + xp2*q1 + xp1*q2 01636 ypp=2*b+xp1*q3+2*xp2*q2+6*dx*q1 01637 01638 return y, yp/dx0, ypp/dx0/dx0 01639 if __name__=="__main__":

---