数值分析作业

最小二乘法

参考：Python闲谈（二）聊聊最小二乘法以及leastsq函数 - 悠望南山 - 博客园 (cnblogs.com)

题目

import matplotlib.pyplot as plt
import numpy as np
import math as math
from scipy.optimize import leastsq

###采样点(Xi,Yi)###
Xi=np.array([0,5,10,15,20,25,30,35,40,45,50,55])
Yi=np.array([0.0,1.27,2.16,2.86,3.44,3.87,4.15,4.37,4.51,4.58,4.62,4.64])

###需要拟合的函数func及误差error###
def func(p,x):
    c,b,a=p
    return c*x**2+b*x+a

def error(p,x,y):
    return func(p,x)-y 

#TEST
p0=[5,2,10]
###主函数从此开始###
Para=leastsq(error,p0,args=(Xi,Yi))
c,b,a=Para[0]
print("c=",c,"b=",b,"a=",a)

plt.scatter(Xi,Yi,color="red",label="Sample Point",linewidth=3) 
x=np.linspace(0,60,1000)
y=c*x**2+b*x+a
plt.plot(x,y,color="orange",label="Fitting Curve",linewidth=2)
plt.legend()
plt.show()

sum=0
for i in range(12):
    sum+=math.pow(error([c,b,a],Xi[i],Yi[i]),2)
print("均方误差为%f" %(sum/12))

题目

import matplotlib.pyplot as plt
import numpy as np
import math as math
from scipy.optimize import leastsq

###采样点(Xi,Yi)###
Xi=np.array([19,25,31,38,44])
Yi=np.array([19.0,32.3,49.0,73.3,97.8])

###需要拟合的函数func及误差error###
def func(p,x):
    b,c,a=p
    return b*x**2+a

def error(p,x,y):
    
    return func(p,x)-y #x、y都是列表，故返回值也是个列表

#TEST
p0=[5,2,10]


###主函数从此开始###
Para=leastsq(error,p0,args=(Xi,Yi)) #把error函数中除了p以外的参数打包到args中
b,c,a=Para[0]
print("b=",b,"a=",a)

plt.scatter(Xi,Yi,color="red",label="Sample Point",linewidth=3) #画样本点
x=np.linspace(0,50,1000)
y=b*x**2+a
plt.plot(x,y,color="orange",label="Fitting Curve",linewidth=2) #画拟合曲线
plt.legend()
plt.show()

sum=0
for i in range(5):
    sum+=math.pow(error([b,0,a],Xi[i],Yi[i]),2)
print("均方误差为%f" %(sum/5))

牛顿插值法

参考：(30条消息) 函数插值法之牛顿插值法 python代码实现_一只懒猫的博客-CSDN博客_python插值函数

题目

import numpy as np
import matplotlib.pyplot as plt

def csb(x,f,j):  
    f0=np.zeros((j+1,x.shape[0]))
    if type(f) is np.ndarray:
        f0[0]=f.copy()
    else:
        for i in range(x.shape[0]):
            f0[0,i]=f(x[i])
    for i in range(1,j+1):
        for k in range(i,j+1):
            f0[i,k]=(f0[i-1,k]-f0[i-1,k-1])/(x[k]-x[k-i])      
    f1=np.vstack([x,f0])
    print('————所求%d阶均差表如下所示————' % j,'\n',f1.T)
    return f1.T

def main():
    x=np.array([0.2,0.4,0.6,0.8,1.0])
    y=np.array([0.98,0.92,0.81,0.64,0.38])
    plt.scatter(x, y, marker='o')
    plt.plot(x, y, label="Newton = 4")
    plt.xlabel("x")
    plt.ylabel("f(x) and Newton(x)")
    plt.ylim(0, 1.20)
    plt.xlim(0, 2.00)
    plt.legend()
    plt.show()
    
    z=csb(x,y,4)
    print("——————预测值——————")
    n=z[0,1]
    X=np.array([0.2,0.28,1.0,1.08])
    for x1 in X:
        n=z[0,1]
        for i in range(4):
            s=1
            k=0
            while k<i+1:
                s=s*(x1-x[k])
                k+=1
            n=n+z[i+1,i+2]*s
        print("P(%f) = %f " % (x1, n))


if __name__ == '__main__':
    main()

龙格函数

参考：拉格朗日插值法 - 简书 (jianshu.com)

题目

from scipy.interpolate import lagrange
import numpy as np
import matplotlib.pyplot as plt

def func(x):
    return 1/(1-25*x**2)

x1=np.linspace(-1,1,10)
y1=func(x1)
lx1 =np.linspace(-1,1,50);
ly1=lagrange(x1,y1)(lx1)

x2=np.linspace(-1,1,20)
y2=func(x2)
lx2 =np.linspace(-1,1,50);
ly2=lagrange(x2,y2)(lx2)

plt.plot(x1, y1, label="n = 10")
plt.plot(x1,y1,"bo", label="point")
plt.plot(lx1, ly1,color="orange", label="lagrange")
plt.legend()
plt.show()
plt.plot(x2, y2, label="n = 20")
plt.plot(x2,y2,"bo", label="point")
plt.plot(lx2, ly2,color="orange", label="lagrange")
plt.legend()
plt.show()

hermite插值

参考：数值计算方法实验之Hermite 多项式插值 (Python 代码) - ぺあ紫泪°冰封ヤ - 博客园 (cnblogs.com)

题目

import numpy as np
from sympy import *
import matplotlib.pyplot as plt


def f(x):
    return  x ** 4


def ff(x):  # f[x0, x1, ..., xk]
    ans = 0
    for i in range(len(x)):
        temp = 1
        for j in range(len(x)):
            if i != j:
                temp *= (x[i] - x[j])
        ans += f(x[i]) / temp
    return ans


def draw(L, newlabel= 'Lagrange插值函数'):
    x = np.linspace(-5, 5, 11)
    y = f(x)
    Ly = []
    for xx in x:
        Ly.append(L.subs(n, xx))
    plt.plot(x, y, "bo",label='point')
    plt.plot(x, Ly, label=newlabel)
    plt.xlabel('x')
    plt.ylabel('y')
    plt.legend()
    plt.show()


def lossCal(L): 
    x = np.linspace(-5, 5, 11)
    y = f(x)
    Ly = []
    for xx in x:
        Ly.append(L.subs(n, xx))
    Ly = np.array(Ly)
    temp = Ly - y
    temp = abs(temp)
    #print(temp.mean())


def calM(P, x):
    Y =   n ** 4
    dfP = diff(P, n)
    return solve(Y.subs(n, x[0]) - dfP.subs(n, x[0]), [m,])[0]


if __name__ == '__main__':
    x = np.array(range(11)) - 5
    y = f(x)
    P = f(x[0])
    
    for i in range(len(x)):
        if i != len(x) - 1:
            temp = ff(x[0:i + 2])
        else:
            temp = m
        for j in x[0:i + 1]:
            temp *= (n - j)
        P += temp
    P = expand(P)
    P = P.subs(m, calM(P, x))
    draw(P, newlabel='Hermite插值多项式')
    lossCal(P)

分段线性插值

参考：[Python] 分段线性插值 - KennyRom - 博客园 (cnblogs.com)

题目

 

import numpy as np
import matplotlib.pyplot as plt
 
#分段线性插值闭包
def get_line(xn, yn):
    def line(x):
        index = -1
         
        #找出x所在的区间
        for i in range(1, len(xn)):
            if x <= xn[i]:
                index = i-1
                break
            else:
                i += 1
         
        if index == -1:
            return -100
         
        #插值
        result = (x-xn[index+1])*yn[index]/float((xn[index]-xn[index+1])) + (x-xn[index])*yn[index+1]/float((xn[index+1]-xn[index]))
         
        return result
    return line
 
xn = [i for i in range(-5,6,1)]
yn = [i**2 for i in xn]
 
#分段线性插值函数
lin = get_line(xn, yn)
 
x = [i for i in range(-5, 6)]
y = [lin(i) for i in x]

#画图
plt.plot(xn ,yn , color="orange",label="original function")
plt.plot(xn, yn, 'ro')
plt.plot(x, y, 'b-',label="segment")
plt.legend()
plt.show()