- 积分
- 129
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2023-3-12
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
# -*- coding: utf-8 -*-
"""
Created on Sat May 6 08:02:43 2023
@author: hp
"""
import numpy as np
import matplotlib.pyplot as plt
#fei0_fei1_fei2
def fei0(x):
return 1
def fei1(x):
return x
def fei2(x):
return pow(x,2)
fig,axes=plt.subplots(figsize=(15,8))
#将函数搞成矩阵形式便于求正规方程组
func=[fei0,fei1,fei2]
x=np.arange(-1,1.2,0.2)
y=1/(1+25*pow(x,2))
#正规方程组
#(fei0,fei0)=fei(0,0)等
fei_fei=np.zeros((3,3))
for i in range(0,3):
for j in range(0,3):
#计算(fei0,fei0)等
for k in range(0,11):
fei_fei[i][j]=fei_fei[i][j]+func[i](x[k])*func[j](x[k])
#计算(fei0,f),(fei1,f)(fei2,f))))
fei_f=np.array([[0.0],[0.0],[0.0]])
for i in range(0,3):
for k in range(0,11):
fei_f[i][0]=fei_f[i][0]+func[i](x[k])*y[k]
a=np.linalg.inv(fei_fei)*fei_f
def nihe(x):
return a[0][0]+a[0][1]*x+a[0][2]*x*x
x2=np.arange(-1,1,0.01)
y2=[0]*len(x2)
y2=nihe(x2)
plt.plot(x2,y2,label="拟合函数")
plt.rcParams['font.sans-serif'] = ['SimHei'] #解决中文显示
plt.rcParams['axes.unicode_minus'] = False #解决符号无法显示
#原函数点
plt.scatter(x,y,c='r',linewidths=5,label="原函数点")
#原函数
y3=1/(1+25*pow(x2,2))
plt.plot(x2,y3,label="原函数",ls='--')
#Newton插值
def newton(x0,y0,x):
#计算差商表并选取其对角线,len(x0)行,len(x0)列
#第0列和第1列已知并输入
a=np.zeros((len(y0)+1,len(x0)))
a[0][0:len(x0)]=np.array(x0)
a[1][0:len(x0)]=np.array(y0)
a=a.T
#先算每一列
for j in range(2,len(y0)+1):
#再算每一行
for i in range(j-1,len(x0)):
a[i][j]=(a[i][j-1]-a[i-1][j-1])/(a[i][0]-a[i+1-j][0])
sum=a[0][1]
#外层循环负责加
for k in range(1,len(x0)):
#内层循环负责实现连乘
t=a[k][k+1]*(x-x0[0])
#从x[1]开始加
for m in range(1,k):
t=t*(x-x0[m])
sum=sum+t
return sum
#x0,y0'
x111=np.arange(-1,1,0.01)
y111=newton(x,y,x111)
plt.plot(x111,y111,label="牛顿插值",ls='-.')
plt.legend(loc=9,fontsize=15)
plt.show()
|
|