• 朝韩将军级会谈时隔11年后在板门店重启 2019-07-23
  • Conférence de presse du Premier ministre chinois 2019-07-09
  • 大观园举办“古琴雅集” 市民感受“古韵端午” 2019-06-30
  • 文脉颂中华——黄河新闻网 2019-06-30
  • 水费欠账竟“穿越”16年?用户质疑:为何没见催缴过? 2019-06-21
  • 买房怎么看风水这个真的实在是太重要了 ——凤凰网房产北京 2019-06-11
  • 6月14日凤凰直通车:茅台再开市场化招聘大门,32个部门要285人葡萄 种植 2019-06-07
  • 习近平为传统文化“代言” 2019-05-27
  • 中巴建交一周年 一系列庆祝活动在巴拿马举行 2019-05-24
  • 招聘启事丨西部网诚聘新媒体编辑记者、实习编辑等人员 2019-05-23
  • A title= href=httpwww.snrtv.comlivech=8 target= 2019-05-23
  • 不止消灭刘海屏 vivo NEX发布会看点汇总 2019-05-22
  • 世相【镜头中的陕西人】 2019-05-20
  • 邓紫棋首任明星制作人 吴亦凡身兼二职 2019-05-20
  • 陶昕然女儿正面照曝光 吃蛋糕萌到爆 2019-05-19
  • Start notebook with --profile=sympy flag.

    Following the 1968 paper by Spath. Consider special case with equally spaced abscissae

    In [1]:
    %qtconsole
    
    In [2]:
    (A, B, C, D)=symbols('A:D') # coefficients
    p = symbols('p') # tension parameter
    y1p = symbols('y1p') # derivative a left endpoint
    ynp = symbols('ynp') # derivative a right endpoint
    xcp = symbols('xcp') # x control point
    ycp = symbols('ycp') # y control point
    

    This is the $k^{th}$ piece of the $n$-piece interpolant,

    $$f_k(x) = A_k+B_k (x-x_k) + C_k \exp(p_k (x-x_k)) +D_k \exp(-p_k (x-x_k))$$

    given the derivative of the target function at $x(0)$ and at the other end $x(n-1)$.

    In [3]:
    f = A(k)+B(k)*(x-x(k)) + C(k)*exp(p(k)*(x-x(k))) +D(k)*exp(-p(k)*(x-x(k)))
    

    Sample data to interpolate

    In [4]:
    X =[0,xcp,1]
    Y =[0,ycp,1]
    

    Set up each piece of interpolant

    In [5]:
    c=[f.subs(x(k),X[i]).subs(k,i) for i in range(2)]
    

    Left-side Interpolation conditions

    In [6]:
    cond_i=[(Y[i]-f.subs(x,x(k)).subs(k,i)) for i in range(2)] # conditions for interpolation
    cond_i+= [ Y[2]- c[1].subs(x,X[2])]
    

    Match end-point 1st derivatives from givens

    In [7]:
    cond_end=[ diff(f,x).subs(k,0).subs(x(0),X[0]).subs(x,X[0]) - y1p,
               diff(f,x).subs(k,1).subs(x(1),X[1]).subs(x,X[2]) - ynp,
     ]
    cond_end
    
    Out[7]:
    [-y1p + B(0) + C(0)*p(0) - D(0)*p(0),
     -ynp + B(1) + C(1)*p(1)*exp((-xcp + 1)*p(1)) - D(1)*p(1)*exp(-(-xcp + 1)*p(1))]

    Inner continuity conditions

    In [8]:
    cond_cont=[] # continuity conditions
    for i,j,v in zip(c[:-2],c[2:],range(1,3)):
        cond_cont.append( (i-j).subs(x,v) )
    cond_cont
    
    cond_cont=[(c[0]-c[1]).subs(x,X[1])]
    

    Inner second derivatives must match

    In [9]:
    d2=[i.diff(x,2) for i in c]
    cond2_cont=[] # 2nd derivative continuity conditions
    for i,j,v in zip(d2[:-2],d2[2:],range(3)):
        cond2_cont.append( (i-j).subs(x,v) )
        
    cond2_cont= [diff(c[0]-c[1],x,2).subs(x,X[1])]
    

    Inner first derivatives must match

    In [10]:
    d=[i.diff(x) for i in c]
    cond1_cont=[] # 2nd derivative continuity conditions
    for i,j,v in zip(d[:-2],d[2:],range(3)):
        cond1_cont.append( (i-j).subs(x,v) )
        
    cond1_cont=[diff(c[0]-c[1],x).subs(x,X[1])]
    
    In [11]:
    len(cond_i)+len(cond_cont)+len(cond_end)+len(cond2_cont)+len(cond1_cont)
    
    Out[11]:
    8
    In [12]:
    for i in (cond_i+cond_cont+cond_end+cond2_cont):
        print i.subs(p(k),0)
    
    -A(0) - C(0) - D(0)
    ycp - A(1) - C(1) - D(1)
    -(-xcp + 1)*B(1) - A(1) - C(1)*exp((-xcp + 1)*p(1)) - D(1)*exp(-(-xcp + 1)*p(1)) + 1
    xcp*B(0) + A(0) - A(1) + C(0)*exp(xcp*p(0)) - C(1) + D(0)*exp(-xcp*p(0)) - D(1)
    -y1p + B(0) + C(0)*p(0) - D(0)*p(0)
    -ynp + B(1) + C(1)*p(1)*exp((-xcp + 1)*p(1)) - D(1)*p(1)*exp(-(-xcp + 1)*p(1))
    C(0)*p(0)**2*exp(xcp*p(0)) - C(1)*p(1)**2 + D(0)*p(0)**2*exp(-xcp*p(0)) - D(1)*p(1)**2
    
    In [13]:
    linsys=(cond_i+cond_cont+cond_end+cond2_cont+cond1_cont)
    M=Matrix([[ l.coeff(i) for i in flatten([[A(j),B(j),C(j),D(j)] for j in range(2)]) ] for l in linsys])
    
    In [14]:
    sum([abs(diff(i,x,2)) for i in c]) # curvature metric
    
    Out[14]:
    Abs(C(0)*p(0)**2*exp(x*p(0)) + D(0)*p(0)**2*exp(-x*p(0))) + Abs(C(1)*p(1)**2*exp((x - xcp)*p(1)) + D(1)*p(1)**2*exp(-(x - xcp)*p(1)))
    In [15]:
    print c[0]
    print c[1]
    
    x*B(0) + A(0) + C(0)*exp(x*p(0)) + D(0)*exp(-x*p(0))
    (x - xcp)*B(1) + A(1) + C(1)*exp((x - xcp)*p(1)) + D(1)*exp(-(x - xcp)*p(1))
    
    In [15]:
     
    
  • 朝韩将军级会谈时隔11年后在板门店重启 2019-07-23
  • Conférence de presse du Premier ministre chinois 2019-07-09
  • 大观园举办“古琴雅集” 市民感受“古韵端午” 2019-06-30
  • 文脉颂中华——黄河新闻网 2019-06-30
  • 水费欠账竟“穿越”16年?用户质疑:为何没见催缴过? 2019-06-21
  • 买房怎么看风水这个真的实在是太重要了 ——凤凰网房产北京 2019-06-11
  • 6月14日凤凰直通车:茅台再开市场化招聘大门,32个部门要285人葡萄 种植 2019-06-07
  • 习近平为传统文化“代言” 2019-05-27
  • 中巴建交一周年 一系列庆祝活动在巴拿马举行 2019-05-24
  • 招聘启事丨西部网诚聘新媒体编辑记者、实习编辑等人员 2019-05-23
  • A title= href=httpwww.snrtv.comlivech=8 target= 2019-05-23
  • 不止消灭刘海屏 vivo NEX发布会看点汇总 2019-05-22
  • 世相【镜头中的陕西人】 2019-05-20
  • 邓紫棋首任明星制作人 吴亦凡身兼二职 2019-05-20
  • 陶昕然女儿正面照曝光 吃蛋糕萌到爆 2019-05-19
  • 一起来捉妖cdkey码领取 白小姐152期资料 pc蛋蛋北京28走势图 中国体彩网大乐透走势 极速时时彩做号 河南快3近200期 广西快乐双彩论坛 北京pk计划人工在线全天免费版 王牌战士腾讯官网测试服名单 博彩老头排列三 3月12日雷霆vs火箭 埃及旋转电子游戏 赫尔南德斯拉齐奥 火影忍者官网 沙特伊蒂哈德鼎盛时期 水果拉霸援彩金