【计算方法】插值和拟合

0.WARNINGS

    本文章主要适用于计算方法代码的实现参考，由于本人是python究极小白，为了实验课速成了一些内容，因此会包含较多的暴力解法orz，如有代码错误、可优化的地方欢迎各位大佬指出，感激不尽。

    此外，本人制作本文的目的，主要是懒的将代码保存在本地中，且博客方便个人的复习。（写这篇的时候已经知道下周要小测了，哇的一声哭出来！）

    插值主要包含的内容有：Lagrange插值（1次线性、2次抛物、n次）；Newton插值（差商）；Hermite插值；分段低次多项式插值&样条插值（高次插值的Runge现象）。

    拟合主要包含的内容有：最小二乘拟合（线性&非线性）。

1.二次Lagrange插值多项式插值法

x0=float(0)
x1=float(1)
x2=float(2)
y0=pow(3,x0)
y1=pow(3,x1)
y2=pow(3,x2)
x=float(input())
f0=(x-x1)*(x-x2)/((x0-x1)*(x0-x2))
f1=(x-x0)*(x-x2)/((x1-x0)*(x1-x2))
f2=(x-x0)*(x-x1)/((x2-x0)*(x2-x1))
ret=f0*y0+f1*y1+f2*y2
print("%.5f" % ret)

2.采用不同阶拉格朗日插值多项式估值

import math

def lagrange(xx,x,n):
    ret=0
    y=list([])
    f=list([])
    for i in range(n+1):
        temp=xx[i]*xx[i]-1
        y.append(pow(math.e,temp))

    for i in range(n+1):
        x1=1
        x2=1
        for j in range(n+1):
            if i!=j:
                x1*=x-xx[j]
                x2*=xx[i]-xx[j]
        f.append(x1/x2)

    for i in range(n+1):
        ret=ret+y[i]*f[i]

    return ret


xx1=input().split()
x=float(input())
xx=list([])
y1=list([])
n=0
for i in xx1:
    xx.append(float(i))
    n=n+1

y1.append(lagrange(xx,x,1))
y1.append(lagrange(xx,x,2))
y1.append(lagrange(xx,x,3))
y1.append(lagrange(xx,x,4))

print("%.4f" % y1[0],' ',"%.4f" % y1[1],' ',"%.4f" % y1[2],' ',"%.4f" % y1[3])

3.用二次插值多项式插值法求f(x)的近似值

import math
x0=float(math.pi/6)
x1=float(math.pi/4)
x2=float(math.pi/3)
y0=0.5
y1=float(math.sqrt(2)/2)
y2=float(math.sqrt(3)/2)

x=float(input())

l0=(x-x1)*(x-x2)/((x0-x1)*(x0-x2))
l1=(x-x0)*(x-x2)/((x1-x0)*(x1-x2))
l2=(x-x1)*(x-x0)/((x2-x1)*(x2-x0))

ret=l0*y0+l1*y1+l2*y2
print("%.4f" % ret)

 4.对于给定的n+1个插值点，计算1到n阶差商值

x0,x1,x2,x3=map(float,input().split())
y0,y1,y2,y3=map(float,input().split())

def f(X1,X2,Y1,Y2):
    return (Y1-Y2)/(X1-X2)

F01=f(x0,x1,y0,y1)
F12=f(x1,x2,y1,y2)
F23=f(x2,x3,y2,y3)

F02=(F01-F12)/(x0-x2)
F13=(F12-F23)/(x1-x3)

F03=(F02-F13)/(x0-x3)

print("%.6f" %(F01))
print("%.6f" %(F02))
print("%.6f" %(F03))
# 暴力解法
# 优化可以采用数组迭代求

 5.用Newton插值公式求f(0.596)的近似值

# 暴力解法，orz
x0=0.40
x1=0.55
x2=0.65
x3=0.80
x4=0.90

y0=0.41075
y1=0.57815
y2=0.69675
y3=0.88811
y4=1.02652

def f(X1,X2,Y1,Y2):
    return (Y1-Y2)/(X1-X2)

F01 = f(x0,x1,y0,y1)
F12 = f(x1,x2,y1,y2)
F23 = f(x2,x3,y2,y3)
F34 = f(x3,x4,y3,y4)

F02 = (F01-F12)/(x0-x2)
F13 = (F12-F23)/(x1-x3)
F24 = (F23-F34)/(x2-x4)

F03 = (F02-F13)/(x0-x3)
F14 = (F13-F24)/(x1-x4)

F04 = (F03-F14)/(x0-x4)

def Newton(x):
    y = y0 + F01*(x-x0) + F02*(x-x0)*(x-x1) + F03*(x-x0)*(x-x1)*(x-x2) + F04*(x-x0)*(x-x1)*(x-x2)*(x-x3)
    return y

X = float(input())
Y = Newton(X)
print("%.3f" %Y)

6.用三次Hermit多项式求 ln(x)的近似值

x0=float(1)
x1=float(2)
y0=float(0)
y1=float(0.693147)
y00=float(1)
y11=float(0.5)

def Hermite(x0,x1,y0,y1,y00,y11,x):
    L0=(x-x1)/(x0-x1)
    L1=(x-x0)/(x1-x0)
    F0=((1+2*L1)*y0+(x-x0)*y00)
    F1=((1+2*L0)*y1+(x-x1)*y11)
    ret=F0*(L0**2)+F1*(L1**2)
    return ret

X=float(input())
ret=Hermite(x0,x1,y0,y1,y00,y11,X)
print("%.5f" %ret)

 7.用线性拟合法求f（x）的近似值

import numpy as np

a0=np.array([[1],[1],[1],[1],[1]])
# fai0
a1=np.array([[0.4],[1.0],[2.0],[4.0],[10.0]])
# fai 1

y=np.array([[0.4053],[1.0071],[2.0167],[3.9963],[10.0165]])
#y

#fai0 * fai0
#fai0 * fai1
#fai1 * fai1
#y * fai0
#y * fai1

b00=0
b01=0
b11=0
r0=0
r1=0

for i in range(0,5):
    b00 += a0[i]*a0[i]
    b01 += a0[i]*a1[i]
    b11 += a1[i]*a1[i]
    r0 += a0[i]*y[i]
    r1 += a1[i]*y[i]

b=(r1*b00-r0*b01) / (b00*b11-b01*b01)
a=(r0-b01*b) / b00

x=float(input())
ret=b*x+a
print("%.2f" % ret)

8.二次曲线拟合数据

import numpy as np

def nonlinear(xi,yi,n):
    a00 = n
    a01 = a02 = a12 = a22 = 0
    for i in range(n):
        a01 += xi[i]
        a02 += xi[i]*xi[i]
        a12 += xi[i]*xi[i]*xi[i]
        a22 += xi[i]*xi[i]*xi[i]*xi[i]
    a11 = a20 = a02
    a21 = a12
    a10 = a01
    a = [[a00, a01,a02], [a10, a11,a12],[a20,a21,a22]]
    a = np.array(a)

    b0 = b1 = b2 = 0
    for i in range(n):
        b0 += yi[i]
        b1 += yi[i] * xi[i]
        b2 += yi[i]*xi[i]*xi[i]
    b = [[b0], [b1], [b2]]
    b = np.array(b)

    x = np.linalg.solve(a, b)
    return x

xi=[1,2,3,4,5,6]
yi=[9.1,18.3,32,47,69.5,94.8]
x0=float(input())
x=nonlinear(xi,yi,6)
y0=x[0]+x[1]*x0+x[2]*x0*x0
print("%.4f" % y0)