最近學(xué)了高等數(shù)值分析，需要做一下數(shù)值分析相關(guān)的編程。感覺(jué)三次樣條插值和Romberg外推加速公式寫(xiě)起來(lái)還是有點(diǎn)難度的。分享一下自己的結(jié)果。

1.三次樣條插值

本來(lái)沒(méi)有什么頭緒，受一個(gè)博主的啟發(fā)，學(xué)習(xí)了他的代碼稍作修改。

原博鏈接：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為二維數(shù)組，用來(lái)存放參數(shù)，size_of_Interval是用來(lái)存放區(qū)間的個(gè)數(shù)
        x_new = np.array(x_vec)
        parameter = []
        size_of_Interval = len(x_new) - 1;
        i = 1
        # 首先輸入方程兩邊相鄰節(jié)點(diǎn)處函數(shù)值相等的方程為2n-2個(gè)方程
        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
        # 輸入端點(diǎn)處的函數(shù)值。為兩個(gè)方程, 加上前面的2n - 2個(gè)方程，一共2n個(gè)方程
        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)
        # 端點(diǎn)函數(shù)一階導(dǎo)數(shù)值相等為n-1個(gè)方程。加上前面的方程為3n-1個(gè)方程。
        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
        # 端點(diǎn)函數(shù)二階導(dǎo)數(shù)值相等為n-1個(gè)方程。加上前面的方程為4n-2個(gè)方程。
        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
        #端點(diǎn)處的函數(shù)值的二階導(dǎo)數(shù)為原函數(shù)的二階導(dǎo)數(shù)，為兩個(gè)方程。總共為4n個(gè)方程。
        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


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

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)


    # 功能：根據(jù)所給參數(shù)，計(jì)算三次函數(shù)的函數(shù)值：
    # 參數(shù):parameters為二次函數(shù)的系數(shù)，x為自變量
    # 返回值：為函數(shù)的因變量
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('三次樣條插值與原函數(shù)的對(duì)比圖', x_axis4, y_axis4)

            
          

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

代碼運(yùn)行結(jié)果

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

            
              import numpy as np

# 編寫(xiě)Romberg求積法，并計(jì)算
def romberg(inf, sup, lenth): #定義函數(shù)的輸入，積分上下界，分塊的數(shù)量
    vec_init = np.zeros(lenth + 1)
    vec_init[0] = 0.5 * (func(sup) + func(inf)) / (sup - inf)
    #初始化T0向量，計(jì)算并賦值
    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
    #設(shè)置停止條件，前后兩次迭代結(jié)果之差小于10^-10
    while np.abs(vec_init[count] - vec_init[count - 1]) > 10 ** -10:
        print('現(xiàn)在在第', 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('計(jì)算結(jié)果為', romberg(0, 1, 20))

#計(jì)算結(jié)果為 0.27219012135491993
            
          

?編寫(xiě)思路：使用一個(gè)向量?jī)?chǔ)存逐次迭代的結(jié)果，比較倒數(shù)1、2位的結(jié)果之差設(shè)為精度條件