2049.sagews

var('x')
def poisci_zlepek(a,b):    #a,b sta seznama
    n=len(a)               #n je stevilo elementov v a seznamu
    h=[]
    u=[]                   #vsi mozni seznami, ki jih potrebujemo
    bi=[]
    v=[]
    z=[]     
    for i in range (0,n-1):  #brez 0!     #napolnimo te sezname
        h.append(a[i+1]-a[i])
        bi.append(6/(h[i])*(b[i+1]-b[i]))
    for j in range (1,n-1):
        u.append(2*(h[j]+h[j-1]))
        v.append(bi[j]-bi[j-1])

    A=matrix(QQ, n-2, n-2, 0)    #matrika je je tista ki ima po diagonali u-je, nad in pod pa h-je, drugje nicle
    #napolnimo matriko A
    for m in range (0,n-2): #napolni brez 0, vkljucno z n-2
        for k in range (0,n-2):
            if (m==k):
                A[m,k]=u[m]
            else:
                if (m+1==k):
                    A[m,k]=h[k]
                else:
                    if (m==k+1):
                        A[m,k]=h[m]

    #napolnili smo jo; sedaj pa resimo sistem -> dobimo z-je
    B=vector(v)           #vektor B je iz v-jev
    Z1=A.solve_right(B)   #je delna resitev, dodati se moramo z0 in zn ki sta oba 0
    z.append(0)           #to so z-ji
    for s in range (0,n-2):
        z.append(Z1[s])
    z.append(0)      #dobili smo z-je
    #napolnimo se A=>Q, B=>W, C=>E
    Q=[]
    W=[]
    E=[]
    for f in range (0,n-1):
        Q.append(1/(6*h[f])*(z[f+1]-z[f])) 
        W.append(z[f]/2)
        E.append(-(h[f]/6)*(z[f+1])-(h[f]/3)*z[f]+(1/h[f])*(b[f+1]-b[f]))
    
    S=[]                        #to bodo naravni kubicni zlepki, kar kot seznam vseh teh
    for p in range (0,n-1):
        S.append(b[p]+(x-a[p])*(E[p]+(x-a[p])*(W[p]+(x-a[p])*Q[p])))

    print html('<FONT color=Red>Naravni kubicni zlepki so:</FONT>')
    for y in range (0,n-1):
        print S[y]

a=[0,1,2,3]
b=[1,1,0,10]
zlepek=poisci_zlepek(a,b)
print zlepek

Naravni kubicni zlepki so:

-(x^2 - 1)*x + 1 (x - 1)*((x - 1)*(4*x - 7) - 2) + 1 -(x - 2)*(3*(x - 5)*(x - 2) - 4) None

a=[-1,0,1]
b=[5,7,9]
zlepek=poisci_zlepek(a,b)
print zlepek

Naravni kubicni zlepki so:

2*x + 7 2*x + 7 None