最近學了高等數值分析，需要做一下數值分析相關的編程。感覺三次樣條插值和Romberg外推加速公式寫起來還是有點難度的。分享一下自己的結果。

1.三次樣條插值

本來沒有什么頭緒，受一個博主的啟發，學習了他的代碼稍作修改。

原博鏈接：https://blog.csdn.net/a19990412/article/details/80574057

            
              import math
import numpy as np
import matplotlib.pyplot  as plt
from sympy import *
from pylab import mpl

def func(y):
    y = np.float64(y)
    return 1/(1 + y * y)

def draw_pic(words, x, y):
    fig=plt.figure()
    plt.plot(x, y, label='interpolation')
    plt.plot(x, func(x), label='raw')
    plt.legend()
    plt.title(words, FontProperties='SimHei')
    #plt.show()
    plt.savefig(words+'.png')
    plt.close(fig)
    pass

def spline3_Parameters(x_vec):
        # parameter為二維數組，用來存放參數，size_of_Interval是用來存放區間的個數
        x_new = np.array(x_vec)
        parameter = []
        size_of_Interval = len(x_new) - 1;
        i = 1
        # 首先輸入方程兩邊相鄰節點處函數值相等的方程為2n-2個方程
        while i < len(x_new) - 1:
            data = np.zeros(size_of_Interval * 4)
            data[(i - 1) * 4] = x_new[i] * x_new[i] * x_new[i]
            data[(i - 1) * 4 + 1] = x_new[i] * x_new[i]
            data[(i - 1) * 4 + 2] = x_new[i]
            data[(i - 1) * 4 + 3] = 1
            data1 = np.zeros(size_of_Interval * 4)
            data1[i * 4] = x_new[i] * x_new[i] * x_new[i]
            data1[i * 4 + 1] = x_new[i] * x_new[i]
            data1[i * 4 + 2] = x_new[i]
            data1[i * 4 + 3] = 1

            parameter.append(data)
            parameter.append(data1)
            i += 1
        # 輸入端點處的函數值。為兩個方程, 加上前面的2n - 2個方程，一共2n個方程
        data = np.zeros(size_of_Interval * 4)
        data[0] = x_new[0] * x_new[0] * x_new[0]
        data[1] = x_new[0] * x_new[0]
        data[2] = x_new[0]
        data[3] = 1
        parameter.append(data)

        data = np.zeros(size_of_Interval * 4)
        data[(size_of_Interval - 1) * 4] = x_new[-1] * x_new[-1] * x_new[-1]
        data[(size_of_Interval - 1) * 4 + 1] = x_new[-1] * x_new[-1]
        data[(size_of_Interval - 1) * 4 + 2] = x_new[-1]
        data[(size_of_Interval - 1) * 4 + 3] = 1
        parameter.append(data)
        # 端點函數一階導數值相等為n-1個方程。加上前面的方程為3n-1個方程。
        i = 1
        while i < size_of_Interval:
            data = np.zeros(size_of_Interval * 4)
            data[(i - 1) * 4] = 3 * x_new[i] * x_new[i]
            data[(i - 1) * 4 + 1] = 2 * x_new[i]
            data[(i - 1) * 4 + 2] = 1
            data[i * 4] = -3 * x_new[i] * x_new[i]
            data[i * 4 + 1] = -2 * x_new[i]
            data[i * 4 + 2] = -1
            parameter.append(data)
            i += 1
        # 端點函數二階導數值相等為n-1個方程。加上前面的方程為4n-2個方程。
        i = 1
        while i < len(x_new) - 1:
            data = np.zeros(size_of_Interval * 4)
            data[(i - 1) * 4] = 6 * x_new[i]
            data[(i - 1) * 4 + 1] = 2
            data[i * 4] = -6 * x_new[i]
            data[i * 4 + 1] = -2
            parameter.append(data)
            i += 1
        #端點處的函數值的二階導數為原函數的二階導數，為兩個方程。總共為4n個方程。
        data = np.zeros(size_of_Interval * 4)
        data[0] = 6 * x_new[0]
        data[1] = 2
        parameter.append(data)
        data = np.zeros(size_of_Interval * 4)
        data[-4] = 6 * x_new[-1]
        data[-3] = 2
        parameter.append(data)
        #df = pd.DataFrame(parameter)
        #df.to_csv('para.csv')
        return parameter


    # 功能：計算樣條函數的系數。
    # 參數：parametes為方程的系數，y為要插值函數的因變量。
    # 返回值：三次插值函數的系數。

def solution_of_equation(parametes, x):
        size_of_Interval = len(x) - 1;
        result = np.zeros(size_of_Interval * 4)
        i = 1
        while i < size_of_Interval:
            result[(i - 1) * 2] = func(x[i])
            result[(i - 1) * 2 + 1] = func(x[i])
            i += 1
        result[(size_of_Interval - 1) * 2] = func(x[0])
        result[(size_of_Interval - 1) * 2 + 1] = func(x[-1])
        result[-2] = 5/13
        result[-1] = -5 / 13
        a = np.array(spline3_Parameters(x))
        b = np.array(result)
        #print(b)
        return np.linalg.solve(a, b)


    # 功能：根據所給參數，計算三次函數的函數值：
    # 參數:parameters為二次函數的系數，x為自變量
    # 返回值：為函數的因變量
def calculate(paremeters, x):
        result = []
        for data_x in x:
            result.append(
                paremeters[0] * data_x * data_x * data_x + paremeters[1] * data_x * data_x + paremeters[2] * data_x +
                paremeters[3])
        return result

x_init4 = np.arange(-5, 5.1, 1 )
result = solution_of_equation(spline3_Parameters(x_init4), x_init4)
#print(spline3_Parameters(x_init4))
#print(result)
x_axis4 = []
y_axis4 = []
for i in range(10):
    temp = np.arange(-5 + i, -4 + i, 0.01)
    x_axis4 = np.append(x_axis4, temp)
    y_axis4 = np.append(y_axis4, calculate(
        [result[4 * i], result[1 + 4 * i], result[2 + 4 * i], result[3 + 4 * i]], temp))
draw_pic('三次樣條插值與原函數的對比圖', x_axis4, y_axis4)

            
          

原博的代碼針對邊界的二次導數為0，故原博使用的矩陣刪去了兩位。不太具有普遍意義，故做了修改。

代碼運行結果

2.Romberg求積分，外推加速公式

            
              import numpy as np

# 編寫Romberg求積法，并計算
def romberg(inf, sup, lenth): #定義函數的輸入，積分上下界，分塊的數量
    vec_init = np.zeros(lenth + 1)
    vec_init[0] = 0.5 * (func(sup) + func(inf)) / (sup - inf)
    #初始化T0向量，計算并賦值
    for i in range(1, lenth):
        vec_init[i] = 0.5 * vec_init[i - 1] + np.array(
            [func(inf + (sup - inf) * (2 * j + 1)/(2 ** i))
             for j in range(2 ** (i - 1) - 1)], dtype=np.float64).sum() * (sup - inf) / 2 ** i
    count = lenth
    deepth = 1
    #設置停止條件，前后兩次迭代結果之差小于10^-10
    while np.abs(vec_init[count] - vec_init[count - 1]) > 10 ** -10:
        print('現在在第', deepth, '層')
        vec_init[0: count - 1] = np.array([(4 ** deepth * vec_init[k + 1] - vec_init[k]) / (4 ** deepth - 1)
                                          for k in range(count - 1)], dtype=np.float64)
        deepth += 1
        count -= 1
    return vec_init[count]

def func(x):
    return (np.log(1 + x)) / (1 + x ** 2)

print('計算結果為', romberg(0, 1, 20))

#計算結果為 0.27219012135491993
            
          

?編寫思路：使用一個向量儲存逐次迭代的結果，比較倒數1、2位的結果之差設為精度條件