第二单元用python学习微积分（十三）牛顿迭代法、中值定理和重要不等式

本文链接：https://blog.csdn.net/bullseye/article/details/122640356

本文内容来自于学习麻省理工学院公开课：单变量微积分-中值定理及重要不等式-网易公开课

开发环境准备：CSDN

二、中值定理（Mean value theorm MVT）

1、罗尔中值定理

(1) 在闭区间 [a,b] 上连续

(2) 在开区间 (a,b) 内可导

(3) f(a)=f(b)，则至少存在一个 c∈(a,b)，使得 f'(c)=0。

2、拉格朗日中值定理(MVT mean value theorem)

一、牛顿迭代法的适用条件

** $x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}$ **

1、误差分析

$E_1 = \left| x - x_1 \right|$

$E_2 = \left| x - x_2 \right|$

...

老师总结是，每迭代一次，解的精度就高一倍

牛顿迭代法要工作的好，需要 $\left| f'(x) \right|$ 不能太小， $\left| f''(x) \right|$ 不能太大（曲线的弯曲程度不能太大，值是0是最适用的情况），x0的取值要在x取值的附近。（这时，老师举了个例子，比如 x^2=5 ,当你想取正数解，但是x0从负数开始给，怎可能得到你不想要的那个解，还有如果选了一个差非常非常远的x0，迭代本身也会遇到问题，所以尽量给算式一个合理的x0的值）

f'(x) 为啥不能太小？

x = symbols('x')
expr =  x**2 -5
def Newton(expr, x0):
    ret = x0 - expr.subs(x, x0)/ expr.diff().subs(x,x0)
    return ret
x1 = Newton(expr, .001)
x2 = Newton(expr, x1)
x3 = Newton(expr, x2)
print('x1 = '+format(float(x1)))
print('x2 = ' +format( float(x2)))
print('x3 = ' +format(float(x3)))

$f(x) = x^2-5$ 可以看到如果 f'(x) = 0 这时迭代式子中分母为0不适用，从几何上看则是函数的切线是横线，永远不会和x轴有交点，即y=0。而如果f'(x)很小，则得到的x1也会离真实值很远，比如上面程序计算的结果，三次迭代还是离真实值很远。

用牛顿迭代解周期性的函数也会出问题，注意给的值不能超过周期

from sympy import *
import numpy as np 

import matplotlib.pyplot as plt 

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.spines['left'].set_position('zero')
ax.spines['bottom'].set_position('zero')
ax.spines['right'].set_color('none')
ax.spines['top'].set_color('none')
ax.xaxis.set_ticks_position('bottom')
ax.yaxis.set_ticks_position('left')
ax.set_aspect(1 ) 

def DrawXY(xFrom,xTo,steps,expr,color,label,plt):
    yarr = []
    xarr = np.linspace(xFrom ,xTo, steps) 
    for xval in xarr:
        yval = expr.subs(x,xval)
        yarr.append(yval)
    y_nparr = np.array(yarr) 
    plt.plot(xarr, y_nparr, c=color, label=label)    

def TangentLine(exprY,x0Val,xVal):
    diffExpr = diff(exprY)
    x1,y1,xo,yo = symbols('x1 y1 xo yo')
    expr = (y1-yo)/(x1-xo) - diffExpr.subs(x,x0Val)
    eq = expr.subs(xo,x0Val).subs(x1,xVal).subs(yo,exprY.subs(x,x0Val))
    eq1 = Eq(eq,0)
    solveY = solve(eq1)
    return xVal,solveY

def DrawTangentLine(exprY, x0Val,xVal1, xVal2, clr, txt):
    x1,y1 = TangentLine(exprY, x0Val, xVal1)
    x2,y2 = TangentLine(exprY, x0Val, xVal2)
    plt.plot([x1,x2],[y1,y2], color = clr, label=txt)
    
def Newton(expr, x0):
    ret = x0 - expr.subs(x, x0)/ expr.diff().subs(x,x0)
    return ret

x = symbols('x')
expr = (sin(x))**2 -1

DrawXY(-4,4,100,y,'blue','y = (sin(x))**2 - 1',plt)

x1 = Newton(expr, 4)
x2 = Newton(expr, x1)
x3 = Newton(expr, x2)
x4 = Newton(expr, x3)
eqy = Eq(expr, 0)
xval = solve(eqy)
print('x1 = '+format(float(x1)))
print('x2 = ' +format( float(x2)))
print('x3 = ' +format(float(x3)))
print('x4 = ' +format(float(x4)))
print('xval[0] = ' +format(float(xval[0])))
print('xval[1] = ' +format(float(xval[1])))
plt.legend(loc='lower right')
plt.show()

另外老师给了一个例子

这个例子会在使用牛顿迭代的情况下，出现循环，注意检验迭代值以避免。

二、中值定理（Mean value theorm MVT）

1、罗尔中值定理

如果 R 上的函数 f(x) 满足以下条件：

(1) 在闭区间 [a,b] 上连续

(2) 在开区间 (a,b) 内可导

(3) f(a)=f(b)，则至少存在一个 c∈(a,b)，使得 f'(c)=0。

证明：因为函数 f(x) 在闭区间[a,b] 上连续，所以存在最大值与最小值，分别用 M 和 m 表示，分两种情况讨论：

a. 若 M=m，则函数 f(x) 在闭区间 [a,b] 上必为常数，结论显然成立。

b. 若 M>m，则因为 f(a)=f(b) 使得最大值 M 与最小值 m 至少有一个在 (a,b) 内某点c处取得，从而c是f(x)的极值点，又条件 f(x) 在开区间 (a,b) 内可导得，f(x) 在 c 处取得极值，由费马引理(费马引理_百度百科)，可导的极值点一定是驻点，推知：f'(c)=0。

2、拉格朗日中值定理(MVT mean value theorem)

如果你从波士顿到洛杉矶(3000mi)，用时6小时，那么某一时刻你能达到整体的平均速度，求这个平均速度。

3000/6 = 500mi/h

如果 f 在区间(a<x<b)内可微, 同时函数在区间( $a\leq x\leq b$ )连续，那么在开区间(a,b)内至少有一点c(a<c<b)使等式 $\frac{f(b) - f(a)}{b-a} = f'(c)$ 成立

(1) 证明

已知f(x) 在(a,b)上连续，在开区间(a,b)内可导，

构造辅助函数 $g(x) = f(x) - f(a) - \frac{f(b)-f(a)}{b-a}(x-a)$

$g(a) = f(a) - f(a) - \frac{f(b)-f(a)}{b-a}(a-a) = 0$

$g(b) = f(b) - f(a) - \frac{f(b)-f(a)}{b-a}(b-a) = 0$

g(a) = g(b)

由于f(x) 在(a,b)上连续，在开区间(a,b)内可导，g(x)则有同样性质

根据罗尔定理可得必有一点 c∈(a,b)，使得 g'(c)=0。

由此可得 $g'(x) = (f(x) - f(a) - \frac{f(b)-f(a)}{b-a}(x-a))' |_{x=c} = 0$

$f'(x) - \frac{f(b)-f(a)}{b-a}(1-0) |_{x=c} = 0$

$f'(x) - \frac{f(b)-f(a)}{b-a} |_{x=c} = 0$

$f'(c) = \frac{f(b)-f(a)}{b-a} c∈(a,b)$

$f(b) -f(a) = (b-a)f'(c)$

$f(b) = f(a) + (b-a)f'(c)$

(2) 几何证明

x = symbols('x')
expr = 5*sin(x)

DrawXY(-4,4,100,expr,'blue','y = 5*sin(x)',plt)
def CalXByDiff(xFrom,xTo,steps,expr,diffValue):
    yarr = []
    xarr = np.linspace(xFrom ,xTo, steps) 
    for xval in xarr:
        yval = expr.diff().subs(x,xval)
        if abs(yval - diffValue)<0.01:
            return xval
    return -1000

xpos = CalXByDiff (-4,4,1000, expr, (expr.subs(x, 2.1) - expr.subs(x, -2.2))/(2.1+2.2))
if(xpos!=-1000):
    DrawTangentLine(expr, xpos, -4,4,'r','f(x) = (f(b)-f(a))/(b-a)')

y = expr.subs(x, 2.1)+((expr.subs(x, 2.1) - expr.subs(x, -2.2))/(2.1+2.2))*(x-2.1)
DrawXY(-4,4,100,y,'green','y = f(a)+(f(b)-f(a)/(b-a))(x-a)',plt)

    
plt.plot([-2.2,-2.2],[-8,6], color = 'gray', label='a = -2.2')
plt.plot([2.1,2.1],[-8,6], color = 'gray', label='b = 2.1')

plt.legend(loc='lower right')
plt.show()