《数学建模简明教程--基于python》学习笔记-第二章-绘图与解方程组-课后习题解答

先纵览一下思考与练习的题目，可以发现主要是画图和解方程方面的题目。

下面是我给出的解答。

准备工作

首先导入几乎每题都需要用的库，并进行绘图和科学计数法的设置。其他模块在使用时再导入。

import matplotlib.pyplot as plt
import numpy as np
plt.rcParams['font.sans-serif']=['SimHei']  # 在图中显示中文字体
plt.rcParams['axes.unicode_minus'] = False  # 在图中显示负号
np.set_printoptions(precision=3, suppress=True) # 设置浮点数小数位数为3 但科学计数法感觉没法关闭

01 绘制双曲函数图像

题1：

x = np.linspace(-10,10,100)

z = np.cosh(x)
plt.plot(x, z, label='双曲余弦函数')

z = np.sinh(x)
plt.plot(x, z, label='双曲正弦函数')

z = np.exp(x)*0.5
plt.plot(x, z, label='y=(1/2)*exp(x)')

plt.axis('square')
plt.xlim([-10,10])
plt.ylim([-10,10])
plt.grid()
plt.legend()

在不同窗口分别绘制以上三个函数，并观察 。

x = np.linspace(-5,5,100)
y = np.cosh(x)

plt.plot(x, y)
plt.title('双曲余弦函数')
print(y.min())  # 注意 min = 1

1.001275651185361

z = np.sinh(x)

plt.plot(x, z)
plt.title('双曲正弦函数')

z = np.exp(x)*0.5

plt.plot(x, z)

​

02 绘制伽马函数图像

题2：

在绘图前你有必要知道伽马函数的一些性质，不然就画个伽马函数会有种一头雾水的感觉
伽马函数可以当成是阶乘在实数集上的延拓，在实数上的伽马函数定义域大于0。
对于正整数n，具有如下性质：
伽马(n) = (n-1) 的阶乘。

下面介绍下伽马函数的由来：
首先,在1728年，哥德巴赫正在考虑数列插值的问题，通俗的说就是把数列的通项公式定义从整数集合延拓到实数集合，例如数列1,4,9,16…

这可以用通项公式n²自然的表达，即便 n 为实数的时候，这个通项公式也是良好定义的。直观的说也就是可以找到一条平滑的曲线y=x²通过所有的整数点(n,n²),从而可以把定义在整数集上的公式延拓到实数集合。

然后有一天，哥德巴赫开始处理阶乘序列1,2,6,24,120,720,…,我们可以计算2!,3!,是否可以计算2.5!呢？我们把最初的一些(n,n!)的点画在坐标轴上，确实可以看到，容易画出一条通过这些点的平滑曲线。
但是哥德巴赫无法解决阶乘往实数集上延拓的这个问题，于是写信请教尼古拉斯·伯努利和他的弟弟丹尼尔·伯努利，由于欧拉当时和丹尼尔·伯努利在一块，他也因此得知了这个问题。欧拉于1729 年解决了这个问题，由此导致了伽玛函数的诞生，当时欧拉只有22岁。 ——百度百科

from scipy.special import gamma  
# SciPy is an open-source software for mathematics, science, and engineering. 

x = np.arange(1,5,0.1)  # 1.0 ~  5.9 共50个点
gamma_z = gamma(x)

plt.plot(x, gamma_z)
plt.axis('square')  # 使得x、y轴显示的单位长度一样长

03 单个窗口绘制二次函数（k=1,2,…,6）

题3:

def fun1(x,k):
    return k*(x**2) + 2*k
    
plt.figure(figsize=(6, 6))
plt.axis('square')  # 使得x,y的单位长度显示比例一致 放后面会出问题

# 绘制图像
for k in range(1,7):
    x = np.linspace(-10, 10, 100)
    z = fun1(x, k) 
    plt.plot(x, z, label = k)
    plt.ylim([0, 20])
    plt.xlim([-10, 10])

# 关于图像的装饰性函数
plt.grid()   
plt.legend()
plt.title('k=1、2、3、4、5、6')
plt.suptitle('y = k*x^2 + 2*k')

04 根据不同K值绘制子图

题4：

row = 2  # 设置子图行数 
col = 3  # 设置子图列数
plt.figure(figsize=(10, 6.18))

for k in range(1,7):
    x = np.linspace(-10, 10, 100)
    z = fun1(x, k) 
    plt.subplot(row, col, k)
    plt.plot(x, z)
    plt.title("k=" + str(k))
    plt.ylim([0,100])
    plt.grid()
    
plt.suptitle('y = k*x^2 + 2*k')

05 绘制二次曲面

题5：

05-1 绘制单叶双曲面

# 单叶双曲面
def fun2(x, y):        
    z = np.sqrt((x**2/4 + y**2/10 - 1)*8)  
    #  z[np.where(np.isnan(z))] = 0  将nan设为0
    return z 

x = np.linspace(-10, 10, 2000)
x, y = np.meshgrid(x, x)  # 有2000*2000个点
z = fun2(x ,y)

fig = plt.figure(figsize=(8,8))
ax = fig.gca(projection='3d')   # 或者 ax = plt.axes(projection='3d')  axes ：axis的复数  这行代码在jupyter中需要与画图在同一格子才起作用
ax.plot_surface(x, y, z, color='y')
ax.plot_surface(x, y, -z, color='y')  # 手动添加开方可取到的负值
plt.show()  # x,y,z的比例只支持auto 不支持equal等

05-2 绘制椭圆双曲面

# 椭圆抛物面
def fun3(z, y): 
    z = x**2/4 + y**2/6
    return z

x = np.linspace(-10, 10, 2000)
x, y = np.meshgrid(x, x)  
z = fun3(x ,y)

fig = plt.figure(figsize=(8,8))
ax = fig.gca(projection='3d')   # 或者 ax = plt.axes(projection='3d')  axes ：axis的复数  这行代码在jupyter中需要与画图在同一格子才起作用
ax.plot_surface(x, y, z, color='y')
plt.show() # x,y,z的比例只支持auto

06 题目无数据跳过

07 求线性方程组的解

题7:

方程组01 （求唯一解）

首先判断方程组解的情况：
发现 R(A) = R(A拔) 所以方程组有唯一解，于是用下面的两种方法求解。

# R（A）=R（A拔）所以有唯一解 
A = np.array([[4,2,-1], 
            [3,-1,2],
            [11,3,0]])
b = np.array([2,10,8])

x = np.linalg.inv(A) @ b
x

array([-4.053e+15,  1.486e+16,  1.351e+16])

# 或者
x = np.linalg.solve(A, b)
x

array([-4.053e+15,  1.486e+16,  1.351e+16])

方程组02 （求最小范数解）

矩阵求秩发现 R(A) = R(A拔) = 2 < 未知变元个数 = 3
所以方程组有无穷多解，这种情况通常求方程组的最小范数解。

A = np.array([[2,3,1],
    [1,-1,2],
    [3,8,-2],
    [4,-1,9]])
b = np.array([4,-15,13,-6])
np.linalg.pinv(A)*b  # 结果是无解情况下的最小二乘解 也是无穷多解情况下的最小范数解 参考 https://blog.csdn.net/weixin_43490741/article/details/104555313 

array([ -25.， 14.， 12.])

08 求非线性方程组的符号解和数值解

题8：

数值解

from scipy.optimize import fsolve

fx = lambda x: [x[0]**2 - x[1] + x[0] -3,
               x[0] + 3*x[1] -2]  # x是列表输入 返回的是“x：”后面的列表

注意：在定义函数时需要把原方程组的右端常数项移到左边

s = fsolve(fx, [1,1])  # [1,1]函数初值
s

array([1.361, 0.213])

符号解

import sympy as sp  

• 符号运算又称计算机代数，就是用计算机推导数学公式
  — 运算以推理方式进行，因此不受截断误差和累计误差的影响
  — 符号计算的速度较慢

x,y = sp.var('x y')  # 使用var函数定义符号变量

sp.solve([x**2 - y + x -3, x + 3*y -2], [x,y]) # 有两个解！ 但数值解只返回了一个

[(-2/3 + sqrt(37)/3, 8/9 - sqrt(37)/9), (-sqrt(37)/3 - 2/3, sqrt(37)/9 + 8/9)]

09 已知f(x)和g(x)的表达式,求非线性方程组的解

题9：

from scipy.optimize import fsolve

def f(x):
    return (abs(x + 1) - abs(x - 1))/2 + np.sin(x)
def g(x):
    return (abs(x + 3) - abs(x - 3))/2 + np.cos(x)

fx = lambda x: [-2*x[0] + 3*f(x[2]) + 4*g(x[3]) - 1,
               -3*x[1] + 2*f(x[2]) + 6*g(x[3]) - 2,
               -x[2] + f(x[0]) + 3*g(x[1]) -3,
               -5*x[3] + 4*f(x[0]) + 6*g(x[1])]  # x是列表输入 返回的是冒号后面的东西

fsolve(fx, [1,1,1,1])  # 这里的[1,1,1,1]将作为求最优解过程中的的初始值

array([4.383, 3.381, 3.14 , 2.477])

10 求超定（矛盾）非线性方程组的最小二乘解

题10：

from scipy.optimize import least_squares

fs = lambda x:[-2*x[0] + 3*f(x[2]) + 4*g(x[3]) - 1,
               -3*x[1] + 2*f(x[2]) + 6*g(x[3]) - 2,
               -x[2]   + f(x[0])   + 3*g(x[1]) - 3,
               -5*x[3] + 4*f(x[0]) + 6*g(x[1]) -1,
              -x[0] - x[2] + f(x[3]) + g(x[1]) - 2,
              -x[1] + 3*x[3] + 2*f(x[0]) - 10*g(x[2]) - 5]

x = least_squares(fx, [0,0,0,0])
x

 active_mask: array([0., 0., 0., 0.])
        cost: 1.575281272016496e-25
         fun: array([-0., -0.,  0.,  0.])
        grad: array([ 0.,  0., -0., -0.])
         jac: array([[-2.   ,  0.   ,  4.731,  3.609],
       [-0.   , -3.   ,  3.154,  5.414],
       [ 1.563,  2.006, -1.   ,  0.   ],
       [ 6.252,  4.012, -0.   , -5.   ]])
     message: '`gtol` termination condition is satisfied.'
        nfev: 8
        njev: 7
  optimality: 4.536140070440464e-12
      status: 1
     success: True
           x: array([-0.973,  0.338, -0.956,  0.098])

学习困惑

1、关于求解非线性方程组，fsolve函数和least_squares函数有什么不同？
2、求解变元个数与方程组个数相同或不同的非线性方程组，应该分别采用什么方法？