计算方法 - 满足第一类边界条件的三次样条插值

题目:
在这里插入图片描述

import numpy as np


def fun(M, X, H, i, x):
    S = M[i-1]*(X[i]-x)**3/(6*H[i]) \
        + M[i]*(x-X[i-1])**3/(6*H[i]) \
        + (Y[i-1]-M[i-1]*H[i]**2/6)*(X[i]-x)/H[i] \
        + (Y[i]-M[i]*H[i]**2/6)*(x-X[i-1])/H[i]
    return S


def make_ans(A, B, C, F):
    n = len(B)
    beta = [C[0]/B[0]]
    for i in range(1, n-1):
        beta.append(C[i]/(B[i]-A[i]*beta[i-1]))
    Y = [F[0]/B[0]]
    for i in range(1, n):
        Y.append((F[i]-A[i]*Y[i-1])/(B[i]-A[i]*beta[i-1]))
    X = np.zeros(n)
    X[n-1] = Y[n-1]
    for i in range(n-2, -1, -1):
        X[i] = Y[i]-beta[i]*X[i+1]
    return X


if __name__ == '__main__':
    X = np.array(input().split(" "))
    Y = np.array(input().split(" "))
    X = X.astype(float)
    Y = Y.astype(float)
    _y0, _yn = input().split(" ")
    _y0 = float(_y0)
    _yn = float(_yn)
    goal_x = float(input())
    n = len(Y)-1
    H = [0]
    for i in range(1, n+1):
        H.append(X[i] - X[i-1])
    Mu = [0]
    Lambda = [1]
    for i in range(1, n):
        Mu.append(H[i]/(H[i] + H[i+1]))
        Lambda.append(1-Mu[i])
    Mu.append(1)
    Lambda.append(0)
    G = [6*((Y[1]-Y[0])/H[1]-_y0)/H[1]]
    for i in range(1, n):
        G.append(6*((Y[i+1]-Y[i])/H[i+1]-(Y[i]-Y[i-1])/H[i])/(H[i]+H[i+1]))
    G.append(6*(_yn-(Y[n]-Y[n-1])/H[n])/H[n])
    G = np.asmatrix(G).T
    B = []
    for i in range(n+1):
        B.append(2)
    M = make_ans(Mu, B, Lambda, G)
    for i in range(n+1):
        print(format(M[i], '.5f'))
    goal_y = 0
    for i in range(n+1):
        if X[i] >= goal_x or i == n:
            goal_y = fun(M, X, H, i, goal_x)
            break
    print(format(goal_y, '.5f'))

    ######################################################
    '''numpy矩阵求逆'''
    W[0][0] = 2
    W[0][1] = 1
    for i in range(1, n):
        W[i][i] = 2
        W[i][i-1] = Mu[i]
        W[i][i+1] = Lambda[i]
    W[n][n] = 2
    W[n][n-1] = 1
    W = np.asmatrix(W)
    M = np.array(W.I * G).reshape(n+1).astype(float)
    ######################################################
  • 0
    点赞
  • 6
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值