数值分析上机实习题
第2章插值法
1. 已知函数在下列各点的值为
试⽤四次⽜顿插值多项式)(x p 4及三次样条韩式)(S x (⾃然边界条件)对数据进⾏插值。⽤图给出(){}10,11,1,0,08.02.0,,x i =+=i x y i i ,)
(x p 4及)(x S Python 代码
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.font_manager import FontProperties
font_set = FontProperties(fname=r"c:\windows\",size=12) #求⽜顿n 次均差 def qiujuncha(x,f,n): for i in range(1,n): for j in range(4,i-1,-1):
f[j]= (f[j] - f[j-1])/(x[j]-x[j-i]) #根据⽜顿多项式求值 def niudun(x,f,x1): sum = f[0]; tmp = 1;
for i in range(1,5): tmp *= (x1-x[i-1]) sum = sum + f[i]*tmp return sum
#⽜顿插值画图 def drawPic(x,f):
x1 = np.linspace(0.2, 1, 100) plt.plot(x1, niudun(x,f,x1))
plt.title(u"⽜顿四次插值",fontproperties=font_set) plt.xlabel(u"x 轴",fontproperties=font_set) plt.ylabel(u"y 轴", fontproperties=font_set) plt.show() def qiu_h(x,h): n = len(x) -1 for i in range(n): print(i)
h[i] = x[i+1]-x[i]
#⾃然边界条件下的三次样条插值求M
def qiu_m(h,f,o,u,d):
n = len(h)
o[0] = 0
u[n] = 0
d[n] = d[0] = 0
a = []
for i in range(1,n):
u[i] = h[i-1]/(h[i-1]+h[i])
for i in range(1,n):
o[i] = h[i]/(h[i-1]+h[i])
for i in range(1,n-1):
d[i] = 6*(f[i+1]-f[i])/(h[i-1]+h[i])
t = [0 for i in range(5)]
t[0] =2
t[1] = o[0]
a.append(t)
for i in range(1,n):
t = [0 for i in range(5)]
t[i - 1] = u [i + 1]
t[i] = 2
t[i + 1] = o [i + 1]
a.append(t)
t = [0 for i in range(5)]
t[n - 1] = u[n]
t[n] = 2
a.append(t)
tmp = np.linalg.solve(np.array(a),np.array(d))
m = []
for i in range(5):
m.append(tmp[i])
return m
#根据三次条插值函数求值
def yangtiao(x1,m,x,y,h,j):
return
m[j]*(x[j+1]-x1)**3/(6*h[j])+m[j+1]*(x1-x[j])**3/(6*h[j])+(y[j]-m[j]*h[j]**2/6)*(x[j+1]-x1)/h[j] +(y[j+1]-m[j+1]*h[j]**2/6)*(x1-x[j])/h[j] def main():
x = [0.2, 0.4, 0.6, 0.8, 1.0]
y = [0.98, 0.92, 0.81, 0.64, 0.38]
f = y[:]
f1 = y[:]
h = [0.2,0.2,0.2,0.2]
u = [0 for n in range(5)]
d = [0 for n in range(5)]
o = [0 for n in range(5)] qiujuncha(x,f,4) qiujuncha(x,f1,2)
m = qiu_m(h,f1,o,u,d) x1 = np.linspace(0.2, 0.4, 10)
p1= plt.plot(x1, yangtiao(x1,m,x,y,h,0),color='red') x1 = np.linspace(0.4, 0.6, 10)
plt.plot(x1, yangtiao(x1, m, x, y, h, 1),color='red') x1 = np.linspace(0.6, 0.8, 10)
plt.plot(x1, yangtiao(x1, m, x, y, h, 2),color='red') x1 = np.linspace(0.8, 1.0, 10)
plt.plot(x1, yangtiao(x1, m, x, y, h, 3),color='red') x1 = np.linspace(0.2, 1.0, 40)
p2 = plt.plot(x1,niudun(x,f,x1),color='green') plt.xlabel(u"x 轴", fontproperties=font_set) plt.ylabel(u"y 轴",
fontproperties=font_set) plt.title("三次样条插值和⽜顿插值")
plt.legend(labels=[u'三次样条插值',u'⽜顿插值'],prop=font_set,loc="best") plt.show() main()
实验结果
运⾏结果可得插值函数图(如图1-1),4次⽜顿插值函数)(x p 4和三次样条插值函数)(x S 如下:
)6.0(*)4.0(*)2.0(625.0)4.0(*)2.0(*3.098.0)(4-------=x x x x x x x P 98.0)8.0(*)6.0(*)4.0(*)2.0(*20833.0+-----x x x x
]4.0,2.0[),2.0(467.4)4.0(9.4)2.0(167.1)(S 3∈-+-+-=x x x x x
]6.0,4.0[),4.0(113.4)6.0(6467.4)4.0(575.1)6.0(167.1)(S 33∈-+-+----=x x x x x x ]8.0,6.0[),6.0(2.3)8.0(113.4)6.0(575.1)(S 3∈-+-+--=x x x x x
]0.1,8.0[),8.0(9.1)0.1(2.3)(S ∈-+-=x x x x
图1-1三次样条插值和⽜顿插值图
2.在区间[-1,1]上分别取n = 10,20⽤两组等距节点对龙格函数做多项式插值三次样条插值,对每个n值画出插值函数及图形。Python代码如下:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.font_manager import FontProperties
font_set = FontProperties(fname=r"c:\windows\",size=12)
def f(x):
return 1/(1+25*x*x)
#求⽜顿n次均差
def qiujuncha(x,f,n):
for i in range(1,n):
for j in range(n-1,i-1,-1):
f[j]= (f[j] - f[j-1])/(x[j]-x[j-i])
#根据⽜顿多项式求值
def niudun(x,f,x1,n):
sum = f[0];
tmp = 1;
for i in range(1,n):
tmp *= (x1-x[i-1])
sum = sum + f[i]*tmp
return sum
def drawPic(n):
linspace numpyx = np.linspace(-1, 1, n)
plt.plot(x, f(x))
plt.show()
def init_x_y(n,x,y):
tmp = 2/n
for i in range(n+1):
x.insert(i, -1 + tmp * i)
y.insert(i,f(-1 + tmp * i))
#⾃然边界条件下的三次样条插值求M def qiu_m(h,f,o,u,d): n = len(h)
o[0] = 0
u[n] = 0
d[n] = d[0] = 0
a = []
for i in range(1,n):
u[i] = h[i-1]/(h[i-1]+h[i])
for i in range(1,n):
o[i] = h[i]/(h[i-1]+h[i])
for i in range(1,n-1):
d[i] = 6*(f[i+1]-f[i])/(h[i-1]+h[i])
t = [0 for i in range(n+1)]
t[0] =2
t[1] = o[0]
a.append(t)
for i in range(1,n):
t = [0 for i in range(n+1)]
t[i - 1] = u [i + 1]
t[i] = 2
t[i + 1] = o [i + 1]
a.append(t)
t = [0 for i in range(n+1)]
t[n - 1] = u[n]
t[n] = 2
a.append(t)
tmp = np.linalg.solve(np.array(a),np.array(d)) m = []
for i in range(n+1):
m.append(tmp[i])
return m
#根据三次条插值函数求值
def yangtiao(x1,m,x,y,h,j):
return
m[j]*(x[j+1]-x1)**3/(6*h[j])+m[j+1]*(x1-x[j])**3/(6*h[j])+(y[j]-m[j]*h[j]**2/6)*(x[j+1]-x1)/h[j] +(y[j+1]-m[j+1]*h[j]**2/6)*(x1-x[j])/h[j] #样条插值函数画图
def draw_yangtiao( m, x, y, h,n):
for i in range(n-1):
x1 = np.linspace(x[i], x[i+1], 10)
plt.plot(x1, yangtiao(x1, m, x, y, h, i),color='red',label=u"三次样条插值") plt.xlabel(u"x轴", fontproperties=font_set)
plt.ylabel(u"y轴", fontproperties=font_set)
plt.title(u"三次样条插值",fontproperties=font_set)
plt.show();
#样条⽜顿函数画图
def draw_niudun(x, f,n):
x1 = np.linspace(-1, 1, 1000)
plt.plot(x1, niudun(x, f, x1,n+1))
print(niudun(x,f,-1,n))
plt.title(u"⽜顿插值", fontproperties=font_set)
plt.xlabel(u"x轴", fontproperties=font_set)
plt.ylabel(u"y轴", fontproperties=font_set)
plt.show()
def main():
n = 20
x=[]

版权声明:本站内容均来自互联网,仅供演示用,请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198,我们将在24小时内删除。