科学计算题目题解

多项式拟合

# encoding: utf-8
from pylab import *
mpl.rcParams['font.sans-serif'] = ['SimHei']
# 中文设置

import numpy as np
import matplotlib.pyplot as plt
import pylab as pl
import matplotlib as mpl

x=[0,34,67,101,135,202,259,336,404,471]
y=[15.18,21.36,25.72,32.29,34.03,39.45,43.15,43.46,40.83,30.75]

plt.plot(x,y,'-',marker="+")

f=np.polyfit(x,y,5)
x2=np.linspace(0,471,300)
y2=np.polyval(f,x2)
plt.plot(x2,y2,'-',color='r')

plt.title('19377XXX XXX')
plt.xlabel('氮肥量')
plt.ylabel('土豆产量')
print(np.polyval(f,300))
plt.show()

从拟合结果来看,拟合多项式次数在4,5次的时候比较接近.当然每个人对拟合效果的评判标准是不一样的,其他的次数也认为是正确的.

核心代码就一句:

f=np.polyfit(x,y,5)

剩下的就是画图了.

拉格朗日插值

# encoding: utf-8
from pylab import *
plt.rcParams['axes.unicode_minus'] = False
mpl.rcParams['font.sans-serif'] = ['SimHei']
# 中文设置

import numpy as np
import matplotlib.pyplot as plt
import pylab as pl
import matplotlib as mpl

def f(x):
    return x/(x*x+1)

def Lagrange(x,y):
    F=np.poly1d([0])
    for i in range(len(x)):
        G=np.poly1d([1])
        for j in range(len(x)):
            if (i!=j):
                G*=(np.poly1d([1,-x[j]])/(x[i]-x[j]))
        F+=G*y[i]
    return F

plt.subplot(2,2,1)
x=np.linspace(-5,5,300)
y=f(x)
plt.plot(x,y,'-',color='r')

plt.subplot(2,2,2)
x1=np.linspace(-5,5,5)
y1=f(x1)
f1=Lagrange(x1,y1)
Y1=np.polyval(f1,x)
plt.plot(x,Y1,'-',color='b')
plt.scatter(x1,y1)

plt.subplot(2,2,3)
x2=np.linspace(-5,5,7)
y2=f(x2)
f2=Lagrange(x2,y2)
Y2=np.polyval(f2,x)
plt.plot(x,Y2,'-',color='green')
plt.scatter(x2,y2)

plt.subplot(2,2,4)
x3=np.linspace(-5,5,10)
y3=f(x3)
f1=Lagrange(x3,y3)
Y3=np.polyval(f1,x)
plt.plot(x,Y3,'-',color='pink')
plt.scatter(x3,y3)

plt.show()

很多人可能对这个题目有一些误解,题目中要求使用拉格朗日插值法构造出一个函数,然后画出这个函数的函数图像.而不是根据题目中给出的插值公式去计算在某一点的值然后再画图.

poly1d函数可以生成一个多项式类型,可以理解为是一种类似int,str的类型.这个类型支持加减乘的运算.在构造的时候,里面传的参数必须是一个list,list的长度决定了生成的多项式次数,里面的数据从前到后对应多项式从高到低次的系数.核心构造代码如下:

def Lagrange(x,y):
    F=np.poly1d([0])
    for i in range(len(x)):
        G=np.poly1d([1])
        for j in range(len(x)):
            if (i!=j):
                G*=(np.poly1d([1,-x[j]])/(x[i]-x[j]))
        F+=G*y[i]
    return F

传入x,y的输入,然后先生成一个0多项式,用来记录求和的结果.通过一个for循环来计算L_{i}(x).原始的公式很长,这里选择拆开计算.每一次只计算其中的一部分,然后累乘起来.在累乘之前最好准备一个常数为1的多项式记录累乘结果.里面的那重循环就是在计算L_{i}(x).计算完毕后按照题意乘以y_{i}累加起来就得到最后的结果了.

在获得了这个多项式后,还需要另外找300以上的点来画图.因为这个题目中一共有四个图,需要对图片位置进行布局.

plt.subplot(2,2,3)

这句的意思是,将画布分成2 \times 2的,然后从左至右,从上至下自动编好号.最后一个参数为3表示,接下来画的图要分配到3号区域.也就是左下角的位置.这样,在每一次画图之前,指定一下区域,就可以实现将四个图画到一个页面上了.

负号显示成方块的,在代码前加一句

plt.rcParams['axes.unicode_minus'] = False

就可以解决.

插值法的选择

# encoding: utf-8
from pylab import *
plt.rcParams['axes.unicode_minus'] = False
mpl.rcParams['font.sans-serif'] = ['SimHei']
# 中文设置

import numpy as np
import matplotlib.pyplot as plt
import pylab as pl
import matplotlib as mpl
import random
import scipy.interpolate as inter

def f(x):
    return x/(x*x+1)

x=np.linspace(-5,5,300)
y=f(x)
plt.plot(x,y,'-',color='r')

x1=np.linspace(-5,5,10)
y1=f(x1)
plt.scatter(x1,y1)

y2=inter.interp1d(x1,f(x1),kind='nearest')
plt.plot(x,y2(x),'-',color='blue')
y3=inter.interp1d(x1,f(x1),kind='linear')
plt.plot(x,y3(x),'-',color='orange')

sum1=sum2=0
for i in range(300):
    w=uniform(-5,5)
    p=f(w)
    q=y2(w)
    t=y3(w)
    sum1+=(p-q)**2
    sum2+=(p-t)**2
print(sum1,sum2)

plt.show()

很多同学写这道题的时候好像忘了上课是讲过现成函数的(我一开始也忘了,写完看标程才发现),直接调函数的话这道题其实就非常简单了,画个图,算个误差就完事了.注意画图的时候,还是要取300个以上的点,如果拿拟合用的那10个点画图的话,就会出现图像重合的问题.最后按照题意随机取点计算误差(拉格朗日插值可以直接套第二题的代码,这里为了简短就省略了).

误差在算完后除不除以n都是可以的,这个题的目的是让大家感受下不同拟合方法的误差的相对大小关系,在评判的时候只要大小关系是对的就没问题,并不会和标程的结果进行比对来给分.而且从结果上可以看出,拉格朗日插值的误差要远大于另外两个,这一点从第二题的最后一个图也能看出来,这就是明显的龙格现象.