第三单元用python学习微积分（二十二）功、平均值、概率（下）和数值积分(1)

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

本文内容来自于学习麻省理工学院公开课：单变量微积分-数值积分-网易公开课

用python解方程组_星-耀的博客-CSDN博客_python解方程组

做功_百度百科

开发环境准备：CSDN

一、概率的经典问题，投飞标问题，问如果一个人向靶心投飞镖，击中的次数与 (c是常数,假设为1)成正比，

首先老师阐释了这个：热量 = 能量 = 做功

热量的单位有千卡(kcal) 能量的单位有千瓦时( kwh ) 做功的单位有焦耳( j )

一、概率的经典问题，投飞标问题，问如果一个人向靶心投飞镖，击中的次数与 $ce^{-r^2}$ (c是常数,假设为1)成正比，

有多大的概率会射到靶子旁边的人。至于这个计算有没有实际意义呢？这里老师提到，二战时有人研究德国的V-2导弹瞄准伦敦发射，会击中哪里，用的就是这个函数 $ce^{-r^2}$ （结果相近！）。

这里击中的概率与距离靶心的长度 r 相关，类正态分布

from sympy import *
import numpy as np 
from scipy.optimize import fsolve
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 DrawInt(xFrom,xTo,steps,expr,color,plt, label=''):
    if(xFrom < 0 and xTo < 0):
        DrawIntNegative(xFrom,xTo,steps,expr,color,plt, label)
    else:
        if(xFrom > 0 and xTo > 0):
            DrawIntPostive(xFrom,xTo,steps,expr,color,plt, label)
        else:
            DrawIntNegative(xFrom,0,steps,expr,color,plt, label)
            DrawIntPostive(0,xTo,steps,expr,color,plt, label)
    
def DrawIntNegative(xFrom1,xTo1,steps,expr,color,plt, label=''):
    xFrom = 0 - xTo1
    xTo = 0 - xFrom1
    width = (xTo - xFrom)/steps
    xarr = []
    yarr = []
    area = 0
    xprev = xFrom
    yvalAll =  0
    xarr.append(0)
    yarr.append(0)
    for step in range(steps):
        yval = expr.subs(x,xprev)
        area += width * yval
        xarr.append(0-xprev)
        yarr.append(0-area)
        xprev= xprev + width
    plt.plot(xarr, yarr, c=color, label =label)
    
def DrawIntPostive(xFrom,xTo,steps,expr,color,plt, label=''):
    width = (xTo - xFrom)/steps
    xarr = []
    yarr = []
    area = 0
    xprev = xFrom
    yvalAll =  0
    xarr.append(0)
    yarr.append(0)
    for step in range(steps):
        yval = expr.subs(x,xprev)
        area += width * yval
        xarr.append(xprev)
        yarr.append(area)
        xprev= xprev + width
    plt.plot(xarr, yarr, c=color, label =label)
    


def DrawParabola(xFrom,xTo,steps,expr,color,plt, label=''):
    width = (xTo - xFrom)/steps
    xarrRect = []
    yarrRect = []
    area = 0
    xprev = xFrom
    yvalAll =  0
    for step in range(int(steps/2)):
        x1 = xprev
        x2 = xprev+width
        x3 = xprev+width+width
        y1 = expr.subs(x,x1) 
        y2 = expr.subs(x,x2)
        y3 = expr.subs(x,x3)
        
        def solveFuncsInDrawParabola(paramlist):
            a,b,c =paramlist[0],paramlist[1],paramlist[2]
            return [a*x1**2+b*x1+c-y1,
                     a*x2**2+b*x2+c-y2,
                     a*x3**2+b*x3+c-y3]
        s=fsolve(solveFuncsInDrawParabola,[0,0,0])
        a,b,c = s[0],s[1],s[2]
        expr1 = a*x**2+b*x+c
        DrawXY(x1,x3,20,expr1,color,'',plt)
        xarrRect.append(x1)
        xarrRect.append(x1)
        xarrRect.append(x3)
        xarrRect.append(x3)
        yarrRect.append(y1)
        yarrRect.append(0)
        yarrRect.append(0)
        yarrRect.append(y3)
        area += 2*width * (y1+4*y2 + y3)/6
        plt.plot(xarrRect, yarrRect, c=color)    
        xprev = x3
    print('============================')    
    if len(label)!=0:
        print(label)    
    print('============================')
    print('width = ', width)
    print('area = ',area)
    areaFinal = (integrate(expr, (x,xFrom,xTo)))
    print ('area final = ',areaFinal)
    print ('ave final = ', areaFinal / (xTo - xFrom))
        
        

def DrawTrapezoidal(xFrom,xTo,steps,expr,color,plt, label=''):
    width = (xTo - xFrom)/steps
    xarrRect = []
    yarrRect = []
    area = 0
    xprev = xFrom
    yvalAll =  0
    for step in range(steps):
        yval1 =expr.subs(x,xprev) 
        yval2 = expr.subs(x,xprev + width)
        xarrRect.append(xprev)
        xarrRect.append(xprev)
        xarrRect.append(xprev + width)
        xarrRect.append(xprev + width)
        xarrRect.append(xprev)
        yarrRect.append(0)
        yarrRect.append(yval1)
        yarrRect.append(yval2)
        yarrRect.append(0)
        yarrRect.append(0)
        area += width * (yval1+yval2)/2
        plt.plot(xarrRect, yarrRect, c=color)    
        xprev= xprev + width
    print('============================')    
    if len(label)!=0:
        print(label)    
    print('============================')
    print('width = ', width)
    print('area = ',area)
    areaFinal = (integrate(expr, (x,xFrom,xTo)))
    print ('area final = ',areaFinal)
    print ('ave final = ', areaFinal / (xTo - xFrom))
        
def DrawRects(xFrom,xTo,steps,expr,color,plt, label=''):
    width = (xTo - xFrom)/steps
    xarrRect = []
    yarrRect = []
    area = 0
    xprev = xFrom
    yvalAll =  0
    for step in range(steps):
        yval = expr.subs(x,xprev + width)
        xarrRect.append(xprev)
        xarrRect.append(xprev)
        xarrRect.append(xprev + width)
        xarrRect.append(xprev + width)
        xarrRect.append(xprev)
        yarrRect.append(0)
        yarrRect.append(yval)
        yarrRect.append(yval)
        yarrRect.append(0)
        yarrRect.append(0)
        area += width * yval
        plt.plot(xarrRect, yarrRect, c=color)    
        xprev= xprev + width
        yvalAll += yval
    print('============================')    
    if len(label)!=0:
        print(label)    
    print('============================')
    print('width = ', width)
    print('ave = ', yvalAll / steps)
    print('area = ',area)
    areaFinal = (integrate(expr, (x,xFrom,xTo)))
    print ('area final = ',areaFinal)
    print ('ave final = ', areaFinal / (xTo - xFrom))

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

def DrawText(x,y, text, fsize=9, clr = 'black'):
    plt.text(x,y, text, fontsize=fsize, color=clr)

c = 1   
x = symbols('x')
expr = c*((np.e)**(-x**2))

DrawXY(-2,2,100,expr,'blue','c*((np.e)**(-x**2))',plt)
plt.plot([1,1],[0,expr.subs(x, 1)], 'red', label='r1')
plt.plot([1.5,1.5],[0,expr.subs(x, 1.5)], 'green', label='r2')

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

上图蓝色曲线围绕 y 轴旋转一周可以形成击中靶子的概率模型，而我们需要考虑的则是r1到r2之间的部分。

看看顶视图（从上向下看）

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) 
x = symbols('x')

expr = (1 - x**2)**0.5

DrawXY(-1,1,100,expr,'red','r1',plt)

expr = -(1 - x**2)**0.5
DrawXY(-1,1,100,expr,'red','',plt)


expr1 = (1.5*1.5 - x**2)**0.5

DrawXY(-1.5,1.5,100,expr1,'green','r2',plt)

expr1 = -(1.5*1.5 - x**2)**0.5
DrawXY(-1.5,1.5,100,expr1,'green','',plt)

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

想象蓝色曲线绕y轴旋转，左侧前视图，右侧顶视图

使用壳层法计算r1到r2之间的体积： $Part = \int_{r_1}^{r_2} 2\pi r \times (height) dr = \int_{r_1}^{r2} 2\pi r \times (ce^{-r^2}) dr$

这里老师一眼看出答案： $-c\pi e^{-r^2}|_{r_1}^{r_2} = c\pi( e^{-r_1^2} - e^{-r_2^2} )$

$Whole = c\pi( e^{-0^2} - e^{-\infty^2} ) ......0 \leq r<\infty$

$= c\pi$

$P(r_1<r<r_2) = \frac{Part}{Whole} = \frac{ c\pi( e^{-r_1^2} - e^{-r_2^2} )}{c\pi} = e^{-r_1^2} - e^{-r_2^2}$

显然这个c确实被削掉了。

同时，把所有的半径都考虑进去，那击中的概率为： $P(0<r<\infty) = 1$

这时老师添加了新的假设：

这个投掷标箭的人在7岁时，靶子的半径为 a , 他中靶的概率为 0.5

$P(0\leq r<a) = \frac{1}{2}$

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) 
x = symbols('x')

a = 1
expr = (a - x**2)**0.5

DrawXY(-a,a,100,expr,'red','r = a',plt)

expr = -(a - x**2)**0.5
DrawXY(-a,a,100,expr,'red','',plt)


r = symbols('r')
xminExpr = r*cos(-30*2/180*np.pi)

expr2 = -((2*a)**2 - x**2)**0.5
x1Min = float(xminExpr.subs(r,2*a)) 
DrawXY(x1Min,2*a,100,expr2,'gray','area where the people stand',plt)

xminExpr = r*cos(-30*2/180*np.pi)
expr3 = -((3*a)**2 - x**2)**0.5
x2Min=float(xminExpr.subs(r,3*a))
DrawXY(x2Min,3*a,100,expr3,'gray','',plt)

plt.plot([2*a, 3*a],[0,0], 'gray', '')
plt.plot([x1Min, x2Min],[expr2.subs(x,x1Min),expr3.subs(x,x2Min)], 'gray', '')

DrawText(3.2*a, 0, '3 o"clock',12)
DrawText(x1Min, expr3.subs(x,x2Min) - 0.2, '5 o"clock',12)

plt.legend(loc='upper left')
plt.show()

由于上图灰色这段环占整个完整的圆环 $\frac{1}{6}$ (....类似2个小时比12个小时...)，所以有

本问题就是求 $\frac{1}{6}P(2a<r<3a)$

由已知： $P(0\leq r<a) = \frac{1}{2}$

$e^{-0^2} - e^{-a^2} = \frac{1}{2}$

$e^{-a^2} = \frac{1}{2}$

$P(2a<r<3a) = e^{-(2a)^2} - e^{-(3a)^2} = e^{-4a^2} - e^{-9a^2} = (e^{-a^2})^4 - (e^{-a^2})^9 = (\frac{1}{2}) ^4 - (\frac{1}{2})^9 \approx (\frac{1}{2}) ^4 =\frac{1}{16}$