[从头学数学] 第169节 函数的应用

剧情提要：

[机器小伟]在[工程师阿伟]的陪同下进入了结丹初期的修炼，

这次要修炼的目标是[函数的应用]。

正剧开始：

星历2016年04月09日 11:38:05, 银河系厄尔斯星球中华帝国江南行省。

[工程师阿伟]正在和[机器小伟]一起研究[函数的应用]。

<span style="font-size:18px;">function funTask(x) {
return Math.log(x)/Math.log(Math.E)+2*x-6;
}

var r = 20;
config.setSector(1,1,1,1);
config.graphPaper2D(0, 0, r);
config.axis2D(0, 0,180, 1);

//坐标轴设定
var scaleX = 2*r, scaleY = 2*r;
var spaceX = 2, spaceY = 2;
var xS = -10, xE = 10;
var yS = -10, yE = 10;
config.axisSpacing(xS, xE, spaceX, scaleX, 'X');
config.axisSpacing(yS, yE, spaceY, scaleY, 'Y');

var transform = new Transform();
//存放函数图像上的点
var a = [], b = [], c = [], d = [], e = [];

//需要显示的函数说明
var f1 = 'y=lnx+2x-6', f2 = 'y = logx', f3 = 'y=x', f4 = 'y=e^x',
f5 = 'y = ln(x)';
//函数描点
for (var x = xS; x <= xE; x+=0.3) {
if (x != 0) {
a.push([x, funTask(x)]);

}
}

//二分法求函数的零点
//区间的最小值和最大值
var minX = 0, maxX = 4;
var y1 = y2 = 0;
x = minX;
y1 = funTask(x);
x = maxX;
y2 = funTask(x);

//如果在给定区间上存在有零点
if (y1 * y2 < 0) {
var epsilon = 0.000001;
while (Math.abs(y1-y2) > epsilon) {

x = minX;
y1 = funTask(x);
x = maxX;
y2 = funTask(x);

x = (minX+maxX)/2;
y = funTask(x);

if (y * y1 < 0) {
maxX = x;
}
else {
minX = x;
}

}
plot.setFillStyle('blue')
.fillText('零点：x = '+ ((minX+maxX)/2).toFixed(3), -200, -60, 120);
}
</span>

这个题，由于[人叫板老师]的精度不够，所以取得的值并不如小伟的准确。

<span style="font-size:18px;">	if (1) {
var r = 20;
config.setSector(1,1,1,1);
config.graphPaper2D(0, 0, r);
config.axis2D(0, 0,180, 1);

//坐标轴设定
var scaleX = 2*r, scaleY = 2*r;
var spaceX = 2, spaceY = 2;
var xS = -10, xE = 10;
var yS = -10, yE = 10;
config.axisSpacing(xS, xE, spaceX, scaleX, 'X');
config.axisSpacing(yS, yE, spaceY, scaleY, 'Y');

var transform = new Transform();
//存放函数图像上的点
var a = [], b = [], c = [], d = [], e = [];

//需要显示的函数说明
var f1 = 'y=2^x+3x-7', f2 = 'y = logx', f3 = 'y=x', f4 = 'y=e^x',
f5 = 'y = ln(x)';
//函数描点
for (var x = xS; x <= xE; x+=0.3) {
if (x != 0) {
a.push([x, funTask(x)]);

}
}

//二分法求函数的零点
//区间的最小值和最大值
var minX = 0, maxX = 2;
var y1 = y2 = 0;
x = minX;
y1 = funTask(x);
x = maxX;
y2 = funTask(x);

//如果在给定区间上存在有零点
if (y1 * y2 < 0) {
var epsilon = 0.000001;
while (Math.abs(y1-y2) > epsilon) {

x = minX;
y1 = funTask(x);
x = maxX;
y2 = funTask(x);

x = (minX+maxX)/2;
y = funTask(x);

if (y * y1 < 0) {
maxX = x;
}
else {
minX = x;
}

}
plot.setFillStyle('blue')
.fillText('零点：x = '+ ((minX+maxX)/2).toFixed(3), -200, -60, 120);
}

//document.write(a+'<br/>');

//存放临时数组
var tmp = [];

//显示变换
if (a.length > 0) {
a = transform.scale(transform.translate(a, 0, 0), scaleX/spaceX, scaleY/spaceY);
//函数1
tmp = [].concat(a);
shape.pointDraw(tmp, 'red');
tmp = [].concat(a);
shape.multiLineDraw(tmp, 'pink');

plot.setFillStyle('red');
plot.fillText(f1, 100, -90, 200);
}

if (b.length > 0) {
b = transform.scale(transform.translate(b, 0, 0), scaleX/spaceX, scaleY/spaceY);
//函数2
tmp = [].concat(b);
shape.pointDraw(tmp, 'blue');
tmp = [].concat(b);
shape.multiLineDraw(tmp, '#22CCFF');

plot.setFillStyle('blue');
plot.fillText(f2, 100, -120, 200);

}

if (c.length > 0) {
c = transform.scale(transform.translate(c, 0, 0), scaleX/spaceX, scaleY/spaceY);
tmp = [].concat(c);
shape.pointDraw(tmp, 'green');
tmp = [].concat(c);
shape.multiLineDraw(tmp, '#CCFF22');
plot.setFillStyle('green');
plot.fillText(f3, 100, -150, 200);
}

if (d.length > 0) {
d = transform.scale(transform.translate(d, 0, 0), scaleX/spaceX, scaleY/spaceY);
tmp = [].concat(d);
shape.pointDraw(tmp, 'orange');
tmp = [].concat(d);
shape.multiLineDraw(tmp, '#CC8800');
plot.setFillStyle('orange');
plot.fillText(f4, 100, -180, 200);
}

if (e.length > 0) {
e = transform.scale(transform.translate(e, 0, 0), scaleX/spaceX, scaleY/spaceY);
tmp = [].concat(e);
shape.pointDraw(tmp, 'purple');
tmp = [].concat(e);
shape.multiLineDraw(tmp, 'purple');
plot.setFillStyle('purple');
plot.fillText(f5, 100, -60, 200);
}

}

}

function funTask(x) {
return Math.pow(2, x)+3*x-7;
}</span>

<span style="font-size:18px;">function funTask(x) {
return Math.pow(x, 3)+1.1*Math.pow(x, 2)+0.9*x-1.4;
}

function funTask(x) {
return 3-Math.log(x)-x;
}
</span>

小伟决定来好好研究一下这个多项式拟合数据点。

这里有一组数据：

<span style="font-size:18px;">//一组点

if (1) {
var r = 20;
config.setSector(1,1,1,1);
config.graphPaper2D(0, 0, r);
config.axis2D(0, 0,180, 1);

//坐标轴设定
var scaleX = 2*r, scaleY = 2*r;
var spaceX = 2, spaceY = 2;
var xS = -10, xE = 10;
var yS = -10, yE = 10;
config.axisSpacing(xS, xE, spaceX, scaleX, 'X');
config.axisSpacing(yS, yE, spaceY, scaleY, 'Y');

var X=[ 1 ,2  ,3 ,4 ,5 ,6],
Y=[ 2.5 ,3.51 ,4.45 ,5.52 ,6.47 ,7.51];

var array = [];
var size = X.length;

for (var i = 0; i < size; i++) {
array.push([X[i], Y[i]]);
}

var transform = new Transform();
var tmp = [];

array = transform.scale(transform.translate(array, 0, 0), scaleX/spaceX, scaleY/spaceY);
tmp = [].concat(array);
shape.pointDraw(tmp, 'green');
//tmp = [].concat(array);
//shape.multiLineDraw(tmp, '#CCFF22');

}</span>

经过这样的拟合计算：

<span style="font-size:18px;">#一次多项式拟合
def linefit(x , y):
N = float(len(x))
sx,sy,sxx,syy,sxy=0,0,0,0,0
for i in range(0,int(N)):
sx  += x[i]
sy  += y[i]
sxx += x[i]*x[i]
syy += y[i]*y[i]
sxy += x[i]*y[i]
a = (sy*sx/N -sxy)/( sx*sx/N -sxx)
b = (sy - a*sx)/N
r = abs(sy*sx/N-sxy)/math.sqrt((sxx-sx*sx/N)*(syy-sy*sy/N))
return a,b,r

def tmp():
X=[ 1 ,2  ,3 ,4 ,5 ,6]
Y=[ 2.5 ,3.51 ,4.45 ,5.52 ,6.47 ,7.51]
a,b,r=linefit(X,Y)
print("X=",X)
print("Y=",Y)
print("拟合结果: y = %10.5f x + %10.5f , r=%10.5f" % (a,b,r) )

>>>
X= [1, 2, 3, 4, 5, 6]
Y= [2.5, 3.51, 4.45, 5.52, 6.47, 7.51]
拟合结果: y =    1.00000 x +    1.49333 , r=   0.99989</span>

来看一下拟合结果是否正确：

<span style="font-size:18px;">	if (1) {
var r = 20;
config.setSector(1,1,1,1);
config.graphPaper2D(0, 0, r);
config.axis2D(0, 0,180, 1);

//坐标轴设定
var scaleX = 2*r, scaleY = 2*r;
var spaceX = 2, spaceY = 2;
var xS = -10, xE = 10;
var yS = -10, yE = 10;
config.axisSpacing(xS, xE, spaceX, scaleX, 'X');
config.axisSpacing(yS, yE, spaceY, scaleY, 'Y');

var X=[ 1 ,2  ,3 ,4 ,5 ,6],
Y=[ 2.5 ,3.51 ,4.45 ,5.52 ,6.47 ,7.51];

var array = [];
var size = X.length;

for (var i = 0; i < size; i++) {
array.push([X[i], Y[i]]);
}

var transform = new Transform();
var tmp = [];

//绘制拟合曲线
var a = [], b = [], c = [], d = [], e = [];

//需要显示的函数说明
var f1 = 'y=x+1.49333';
//函数描点
for (var x = xS; x <= xE; x+=1) {
if (x != 0) {
a.push([x, funTask(x)]);
}
}

if (a.length > 0) {
a = transform.scale(transform.translate(a, 0, 0), scaleX/spaceX, scaleY/spaceY);
//函数1
tmp = [].concat(a);
shape.pointDraw(tmp, 'red');
tmp = [].concat(a);
shape.multiLineDraw(tmp, 'pink');

plot.setFillStyle('red');
plot.fillText(f1, 100, -90, 200);
}

//原数据点
array = transform.scale(transform.translate(array, 0, 0), scaleX/spaceX, scaleY/spaceY);
tmp = [].concat(array);
shape.pointDraw(tmp, 'green');
//tmp = [].concat(array);
//shape.multiLineDraw(tmp, '#CCFF22');

}

function funTask(x) {
return x+1.49333;
}</span>

好像是没错。

来看看二次多项式拟合。

<span style="font-size:18px;">def polyfit(x, y, degree):
import numpy;
results = {}
coeffs = numpy.polyfit(x, y, degree)
results['polynomial'] = coeffs.tolist()

# r-squared
p = numpy.poly1d(coeffs)
# fit values, and mean
yhat = p(x)                         # or [p(z) for z in x]
ybar = numpy.sum(y)/len(y)          # or sum(y)/len(y)
ssreg = numpy.sum((yhat-ybar)**2)   # or sum([ (yihat - ybar)**2 for yihat in yhat])
sstot = numpy.sum((y - ybar)**2)    # or sum([ (yi - ybar)**2 for yi in y])
results['determination'] = ssreg / sstot #准确率
return results

def tmp():
x=[ 1 ,2  ,3 ,4 ,5 ,6]
y=[ 2.5 ,3.51 ,4.45 ,5.52 ,6.47 ,7.2]
z1 = polyfit(x, y, 2)
print(z1);

>>>
{'determination': 0.99893014166855909, 'polynomial': [-0.024285714285714542, 1.1257142857142872, 1.3699999999999986]}</span>

<span style="font-size:18px;">	if (1) {
var r = 20;
config.setSector(1,1,1,1);
config.graphPaper2D(0, 0, r);
config.axis2D(0, 0,180, 1);

//坐标轴设定
var scaleX = 2*r, scaleY = 2*r;
var spaceX = 2, spaceY = 2;
var xS = -10, xE = 10;
var yS = -10, yE = 10;
config.axisSpacing(xS, xE, spaceX, scaleX, 'X');
config.axisSpacing(yS, yE, spaceY, scaleY, 'Y');

var X=[ 1 ,2  ,3 ,4 ,5 ,6]
Y=[ 2.5 ,3.51 ,4.45 ,5.52 ,6.47 ,7.2]

var array = [];
var size = X.length;

for (var i = 0; i < size; i++) {
array.push([X[i], Y[i]]);
}

var transform = new Transform();
var tmp = [];

//绘制拟合曲线
var a = [], b = [], c = [], d = [], e = [];

//需要显示的函数说明
var f1 = 'y=-0.024x2+1.126x+1.370';
//函数描点
for (var x = xS; x <= xE; x+=1) {
if (x != 0) {
a.push([x, funTask(x)]);
}
}

if (a.length > 0) {
a = transform.scale(transform.translate(a, 0, 0), scaleX/spaceX, scaleY/spaceY);
//函数1
tmp = [].concat(a);
shape.pointDraw(tmp, 'red');
tmp = [].concat(a);
shape.multiLineDraw(tmp, 'pink');

plot.setFillStyle('red');
plot.fillText(f1, 100, -90, 200);
}

//原数据点
array = transform.scale(transform.translate(array, 0, 0), scaleX/spaceX, scaleY/spaceY);
tmp = [].concat(array);
shape.pointDraw(tmp, 'green');
//tmp = [].concat(array);
//shape.multiLineDraw(tmp, '#CCFF22');

}
</span>

这也不知道它怎么算的，明显有很大偏差啊。

这里顺便贴下对数拟合的代码，这个和小伟没任何关系，

因为不是小伟写的，一时也不好验证对错，先放着看看吧。

<span style="font-size:18px;">对数函数拟合-这个是最难的，baidu上都找不到，google了半天才找到的。指数、幂数拟合啥的，都用这个，把func改写一下就行
from scipy import log as log print pcov
import numpy
from scipy import log
from scipy.optimize import curve_fit

def func(x, a, b):
y = a * log(x) + b
return y

def polyfit(x, y, degree):
results = {}
#coeffs = numpy.polyfit(x, y, degree)
popt, pcov = curve_fit(func, x, y)
results['polynomial'] = popt

# r-squared
yhat = func(x ,popt[0] ,popt[1] )                         # or [p(z) for z in x]
ybar = numpy.sum(y)/len(y)          # or sum(y)/len(y)
ssreg = numpy.sum((yhat-ybar)**2)   # or sum([ (yihat - ybar)**2 for yihat in yhat])
sstot = numpy.sum((y - ybar)**2)    # or sum([ (yi - ybar)**2 for yi in y])
results['determination'] = ssreg / sstot

return results

x=[ 1 ,2  ,3 ,4 ,5 ,6]
y=[ 2.5 ,3.51 ,4.45 ,5.52 ,6.47 ,7.51]
z1 = polyfit(x, y, 2)
print z1
</span>

本节到此结束，欲知后事如何，请看下回分解。