第二章例题与课后作业
例2.1字符串操作示例
点击查看代码
str1 = "Hello World!"
print(str1)
print(str1[0:-1])
print(str1[-1])
print(str1[2:5])
print(str1[2:])
print(str1 * 2)
print("学号:3001")
例2.2统计下列5行字符串中字符a,c,g,t出现的的频数
点击查看代码
import numpy as np
a = []
with open("D:\数学建模\例2.2.txt") as f:
for (i, s) in enumerate(f):
a.append([s.count('a'), s.count('c'),
s.count('g'),s.count('t')])
b = np.array(a); print(b)
print("学号:3001")
例2.3列表操作示例
点击查看代码
L = ['abc', 12, 3.45, 'Python', 2.789]
print(L)
print(L[0])
L[0] = 'a'
L[1:3] = ['b', 'Hello']
print(L)
L[2:4] = []
print(L)
print("学号:3001")
例2.4使用列表推导式实现嵌套列表的平铺
点击查看代码
a = [[1,2,3],[4,5,6],[7,8,9]]
d = [c for b in a for c in b]
print(d)
print("学号:3001")
例2.5在列表推导式中使用if过滤不符合条件的元素
点击查看代码
import os
fn = [filename for filename in
os.listdir('D:\PythonWork')
if filename.endswith(('.exe','.py'))]
print(fn)
print("学号:3001")
from numpy.random import randint
import numpy as np
a = randint(10,20,16)
ma = max(a)
ind1 = [index for index,value in enumerate(a) if value == ma]
ind2 = np.where(a == ma)
print(ind1);print(ind2[0])
print("学号:3001")
例2.6元组操作示例
点击查看代码
T = ('abc', 12, 3.45, 'Python', 2.789)
print(T)
print(T[-1])
print(T[1:3])
print("学号:3001")
例2.7集合操作示例
点击查看代码
student = {'Tom', 'Jim', 'Mary', 'Tom', 'Jack', 'Rose'}
print(student)
a = set('abcdabc')
print(a)
print("学号:3001")
例2.8字典操作示例
点击查看代码
dict1 ={'Alice': '123', 'Beth': '456', 'Cecil': 'abc'}
print(dict1['Alice'])
dict1['new'] = 'Hello'
dict1['Alice'] = '1234'
dict2 = {'abc': 123, 456: 78.9}
print(dict2[456])
print("学号:3001")
例2.9字典的get()方法使用示例
点击查看代码
Dict={'age':18,'score':[98,97],'name':'Zhang','sex':'male'}
print(Dict['age']) #输出18
print(Dict.get('age')) #输出18
print(Dict.get('address','Not Exists.'))#输出 No Exists.
print(Dict['address']) #出错
print("学号:3001")
例2.10字典元素的访问示例
点击查看代码
Dict={'age':18,'score':[98,97],'name':'Zhang','sex':'male'}
for item in Dict:
print(item)
print("----------")
for item in Dict.items():
print(item)
print("----------")
for value in Dict.values():
print(value)
print("学号:3001")
例2.11首先生成包含1000个随机字符的字符串,然后统计每个字符的出现次数,注意get()方法的应用
2.11的(1)
点击查看代码
import string
import random
x=string.ascii_letters+string.digits
y=''.join([random.choice(x) for i in range(1000)])
#choice()用于从多个元素中随机选择一个
d=dict() #构造空字典
for ch in y:
d[ch]=d.get(ch,0)+1;
for k,v in sorted(d.items()):
print(k,':',v)
print("学号:3001")
2.11的(2)
点击查看代码
import string, random, collections #依次加载三个模块
x=string.ascii_letters+string.digits
y=''.join([random.choice(x) for i in range(1000)])
count=collections.Counter(y)
for k,v in sorted(count.items()):
print(k, ':', v)
print("学号:3001")
例2.12分别编写求n!和输出斐波那契数列的函数,并调用两个函数进行测试
点击查看代码
def factorial(n):
r = 1
while n > 1:
r *= n
n -= 1
return r
def fib(n):
a, b = 1, 1
while a < n:
print(a, end=' ')
a, b = b, a+b
print('%d! = %d'%(5,factorial(5)))
fib(200)
print("学号:3001")
例2.13数据分组
点击查看代码
def bifurcate_by(L, fn):
return [[x for x in L if fn(x)],
[x for x in L if not fn(x)]]
s=bifurcate_by(['beep', 'boop', 'foo', 'bar'], lambda x: x[0] =='b')
print(s)
print("学号:3001")
例2.14用匿名函数,求3个数的乘积和列表元素的值
点击查看代码
f=lambda x, y, z: x * y * z
L=lambda x: [x**2, x**3, x**4]
print(f(3,4,5)); print(L(2))
print("学号:3001")
例2.15加载模块示例
点击查看代码
import math
import random
import numpy.random as nr
a=math.gcd(12,21)
b=random.randint(0,2)
c=nr.randint(0,2,(4,3))
print(a); print(b); print(c)
print("学号:3001")
例2.16导入模块示例
点击查看代码
from random import sample
from numpy.random import randint
a=sample(range(10),5)
b=randint(0,10,15)
print(a); print(b)
print("学号:3001")
例2.17导入模块示例
点击查看代码
from math import *
a=sin(3)
b=pi
c=e
d=radians(180)
print(a);print(b);print(c);print(d)
print("学号:3001")
例2.18调入自定义函数factorial()和fin(),计算6!,输出300以内的斐波那契数列
点击查看代码
def factorial(n):
r = 1
while n > 1:
r *= n
n -= 1
return r
def fib(n):
a,b = 1,1
while a < n:
print(a,end=' ')
a,b = b,a+b
print('%d! = %d'%(6,factorial(6)))
fib(300)
print("学号:3001")
例2.19sorted()使用示例
点击查看代码
import numpy.random as nr
x1 = list(range(9, 21))
nr.shuffle(x1) # shuffle()用来随机打乱顺序
x2 = sorted(x1) # 按照从小到大排序
x3 = sorted(x1, reverse = True) # 按照从大到小排序
x4 = sorted(x1, key = lambda item: len(str(item))) # 以指定的规则排序
print(x1); print(x2); print(x3); print(x4)
print("学号:3001")
例2.20enumerate()函数使用示例
点击查看代码
x1 = 'abcde'
x2 = list(enumerate(x1))
for ind, ch in enumerate(x1):
print(ch)
print("学号:3001")
例2.21map()函数使用示例
点击查看代码
import random
x = random.randint(1e5, 1e8) # 生成一个随机整数
y = list(map(int, str(x))) # 提出每位上的数字
z = list(map(lambda x, y: x%2 == 1 and y%2 == 0, [1, 3, 2, 4, 1], [3, 2, 1, 2]))
print(x); print(y); print(z)
print("学号:3001")
例2.22filter()函数使用示例
点击查看代码
a = filter(lambda x: x>0, [1, 11, 2, 45, 7, 6, 13])
b = filter(lambda x: x.isalnum(), ['abc', 'xy12', '***'])
# isalnum() 是测试是否为字母或数字的方法
print(list(a)); print(list(b))
print("学号:3001")
例2.23过滤重复值
点击查看代码
def filter_non_unique(L):
return [item for item in L if L.count(item) == 1]
a = filter_non_unique([1, 2, 2, 3, 4, 4, 5])
print(a)
print("学号:3001")
例2.24zip()函数使用示例
点击查看代码
s1 = [str(x)+str(y) for x, y in zip(['v']*4, range(1, 5))]
s2 = list(zip('abcd', range(4)))
print(s1); print(s2)
print("学号:3001")
例2.25数组生成示例1
点击查看代码
import numpy as np
a1 = np.array([1, 2, 3, 4]) # 生成整型数组
a2 = a1.astype(float)
a3 = np.array([1, 2, 3, 4], dtype = float) # 浮点数
print(a1.dtype); print(a2.dtype); print(a3.dtype)
b = np.array([[1, 2, 3], [4, 5, 6]])
c = np.arange(1, 4) # 生成数组 [1, 2, 3, 4]
d = np.linspace(1, 4, 4) # 生成数组 [1, 2, 3, 4]
e = np.logspace(1, 3, 3, base = 2) # 生成数组 [2, 4, 8]
print(b);print(c);print(d);print(e)
print("学号:3001")
例2.26数组生成示例2
点击查看代码
import numpy as np
a = np.ones(4, dtype = int) # 输出[1, 1, 1, 1]
b = np.ones((4,), dtype = int) # 同a
c = np.ones((4, 1)) # 输出 4 行 1 列的数组
d = np.zeros(4) # 输出 [0, 0, 0, 0]
e = np.empty(3) # 生成 3 个元素的空数组行向量
f = np.eye(3) # 生成 3 阶单位矩阵
g = np.eye(3, k = 1) # 生成第 k 对角线元素为 1,其它元素为 0 的 3 阶方针
h = np.zeros_like(a) # 生成与 a 同维的全 0 数组
print(a);print(b);print(c);print(d)
print(e);print(f);print(g);print(h)
print("学号:3001")
例2.27数组元素的索引示例
点击查看代码
import numpy as np
a = np.arange(16).reshape(4,4) # 生成 4 行 4 列的数组
b = a[1][2] # 输出6
c = a[1, 2] # 同 b
d = a[1:2, 2:3] # 输出[[6]]
x = np.array([0, 1, 2, 1])
print(a[x==1]) # 输出 a 的第 2、4 行元素
print("学号:3001")
例2.28矩阵合并示例
点击查看代码
import numpy as np
a = np.arange(16).reshape(4, 4) # 生成 4 行 4 列的数组
b = np.floor(5*np.random.random((2, 4)))
c = np.ceil(6*np.random.random((4, 2)))
d = np.vstack([a, b]) # 上下合并矩阵
e = np.hstack([a, c]) # 左右合并矩阵
print(a);print(b);print(c);print(d);print(e)
print("学号:3001")
例2.29矩阵分割示例
点击查看代码
import numpy as np
a = np.arange(16).reshape(4, 4) # 生成 4 行 4 列的数组
b = np.vsplit(a, 2) # 行分割
print('行分割:\n', b[0], '\n', b[1])
c = np.hsplit(a, 4) # 列分割
print('列分割:\n', c[0], '\n', c[1], '\n', c[2], '\n', c[3])
print("学号:3001")
例2.30矩阵元素求和示例
点击查看代码
import numpy as np
a = np.array([[0, 3, 4], [1, 6, 4]])
b = a.sum() # 使用方法,求矩阵所有元素的和
c1 = sum(a) # 内置函数求矩阵逐列元素的和
c2 = np.sum(a, axis=0) # 使用函数,求矩阵逐列元素的和
c3 = np.sum(a, axis=0, keepdims=True) # 逐列求和
print(c2.shape, c3.shape) # c2是(3,)数组,c3是(1, 3)数组
print("学号:3001")
例2.31逐个元素运算示例
点击查看代码
import numpy as np
a = np.array([[0, 3, 4], [1, 6, 4]])
b = np.array(([1, 2, 3], [2, 1, 4]))
c = a/b # 两个矩阵对应元素相除
d = np.array([2, 3, 2])
e = a*b # d 先广播成与 a 同维数的矩阵,再逐个元素相乘
f = np.array([[3], [2]])
g = a*f # f 先广播成与 a 同维数的矩阵,再逐个元素相乘
h = a**(1/2) # a 矩阵逐个元素的 1/2 次幂
print(a);print(b);print(c);print(d)
print(e);print(f);print(g);print(h)
print("学号:3001")
例2.32矩阵乘法示例
点击查看代码
import numpy as np
a = np.ones(4)
b = np.arange(2, 10, 2)
c = a @ b # a 作为行向量,b 作为列向量
d = np.arange(16).reshape(4, 4)
f = a @ d # a 作为行向量
g = d @ a # a 作为列向量
print(a);print(b);print(c)
print(d);print(f);print(g)
print("学号:3001")
例2.33求下列矩阵的各个行向量的2范数,各个列向量的2范数和矩阵2范数
点击查看代码
import numpy as np
a = np.array([[0, 3, 4], [1, 6, 4]])
b = np.linalg.norm(a, axis = 1) # 求行向量 2 范数
c = np.linalg.norm(a, axis = 0) # 求列向量 2 范数
d = np.linalg.norm(a) # 求矩阵 2 范数
print('行向量 2 范数为:', np.round(b, 4))
print('列向量 2 范数为:', np.round(c, 4))
print('矩阵 2 范数为:', np.round(d, 4))
print("学号:3001")
例2.34求解线性方程组
点击查看代码
import numpy as np
a = np.array([[3, 1],[1, 2]])
b = np.array([9, 8])
x1 = np.linalg.inv(a) @ b # 第一种解法
x2 = np.linalg.solve(a, b) # 第二种解法
print(x1); print(x2)
print("学号:3001")
例2.35求解线性方程组
点击查看代码
import numpy as np
a = np.array([[3, 1], [1, 2], [1, 1]])
b = np.array([9, 8, 6])
x = np.linalg.pinv(a) @ b
print(np.round(x, 4))
print("学号:3001")
例2.36求下列矩阵的特征值和特征向量
点击查看代码
import numpy as np
a = np.eye(4)
b = np.rot90(a)
c, d = np.linalg.eig(b)
print('特征值为:', c)
print('特征向量为:\n', d)
print("学号:3001")
例2.37生成服从标准正态分布的24*4随机数矩阵,并保持为DataFrame数据结构
点击查看代码
import pandas as pd
import numpy as np
dates = pd.date_range(start = '20191101', end = '20191124', freq = 'D')
a1 = pd.DataFrame(np.random.randn(24, 4), index = dates, columns = list('ABCD'))
a2 = pd.DataFrame(np.random.rand(24, 4))
print(a1);print(a2)
print("学号:3001")
例2.38数据写入文件示例
点击查看代码
import pandas as pd
import numpy as np
# 生成日期索引和随机数据
dates = pd.date_range(start='20191101', end='20191124', freq='D')
a1 = pd.DataFrame(np.random.randn(24, 4), index=dates, columns=list('ABCD'))
a2 = pd.DataFrame(np.random.randn(24, 4))
# 将DataFrame保存到Excel和CSV文件
a1.to_excel('data2_38_1.xlsx')
a2.to_csv('data2_38_2.csv')
# 使用with语句自动管理ExcelWriter对象的生命周期
with pd.ExcelWriter('data2_38_3.xlsx') as f:
a1.to_excel(f, "Sheet1")
a2.to_excel(f, "Sheet2")
print("学号:3001")
例2.39从文件中读入数据示例
点击查看代码
import pandas as pd
a = pd.read_csv('data2_38_2.csv', usecols = range(1, 5))
b = pd.read_excel('data2_38_3.xlsx', "Sheet2", usecols = range(1, 5))
print(a); print(b)
print("学号:3001")
例2.40DataFrame数据的拆分、合并和分组计算示例
点击查看代码
import pandas as pd
import numpy as np
d = pd.DataFrame(np.random.randint(1, 6, (10, 4)), columns = list('ABCD'))
d1 = d[:4] # 获取前 4 行数据
d2 = d[4:] # 获取第 5 行以后的数据
dd = pd.concat([d1, d2]) # 数据行合并
s1 = d.groupby('A').mean() # 数据分组求均值
s2 = d.groupby('A').apply(sum) # 数据分组求和
print(d1);print(d2)
print("学号:3001")
例2.41DataFrame数据操作示例
点击查看代码
import pandas as pd
import numpy as np
a = pd.DataFrame(np.random.randint(1, 6, (5, 3)),
index = ['a', 'b', 'c', 'd', 'e'],
columns = ['one', 'two', 'three'])
a.loc['a', 'one'] = np.nan # 修改第 1 行第 1 列的数据
b = a.iloc[1:3, 0:2].values # 提取第 2、3 行,第 1、2 列数据
a['four'] = 'bar' # 增加第 4 列数据
a2 = a.reindex(['a', 'b', 'c', 'd', 'e', 'f'])
a3 = a2.dropna() # 删除有不确定值的行
print(a2);print(a3)
print("学号:3001")
例2.42遍历文件data2_2.txt中所有行,统计每一行中字符的个数
点击查看代码
with open('data2_2.txt',encoding='utf-8') as fp:
L1=[]; L2=[];
for line in fp:
L1.append(len(line))
L2.append(len(line.strip())) # 去掉换行符
data = [str(num)+'\t' for num in L2] # 转换为字符串
print(L1); print(L2)
with open('data2_2.txt', 'w') as fp2:
fp2.writelines(data)
print("学号:3001")
例2.43随机产生一个数据矩阵,把它存入具有不同分隔符格式的文本文件中,再把数据从文本文件中提取出来
点击查看代码
import numpy as np
a = np.random.rand(6, 8) # 生成 6*8 的[0, 1)上均匀分布的随机数矩阵
np.savetxt("data2_43_1.txt", a) # 存在以制表符分隔的文本文件
np.savetxt("data2_43_2.csv", a, delimiter=',') # 存成以逗号分隔的 CSV 文件
b = np.loadtxt("data2_43_1.txt") # 加载空格分隔的文本文件
c = np.loadtxt("data2_43_2.csv", delimiter=',') # 加载 CSV 文件
print("学号:3001")
例2.44求方程在给定初值1.5附近的一个实根
点击查看代码
from scipy.optimize import fsolve, root
fx = lambda x: x**980-5.01*x**979+7.398*x**978\
-3.388*x**977-x**3+5.01*x**2-7.398*x+3.388
x1 = fsolve(fx, 1.5, maxfev=4000) #函数调用4000次
x2 = root(fx, 1.5)
print(x1,'\n','-------------'); print(x2)
print("学号:3001")
例2.45求下列方程组的一组数值解
点击查看代码
from scipy.optimize import fsolve, root
fx = lambda x: [x[0]**2+x[1]**2-1, x[0]-x[1]]
s1 = fsolve(fx, [1, 1])
s2 = root(fx, [1, 1])
print(s1,'\n','--------------'); print(s2)
print("学号:3001")
例2.46分别计算a=2,b=1; a=2,b=10时,I(a,b)的值
点击查看代码
from scipy.integrate import quad
def fun42(x, a, b):
return a*x**2+b*x
I1 = quad(fun42, 0, 1, args=(2, 1))
I2 = quad(fun42, 0, 1, args=(2, 10))
print(I1); print(I2)
print("学号:3001")
例2.47已知4个观测站的位置坐标(xi,yi)(i=1,2,3,4),每个观测站探测到距未知信号的距离di(i=1,2,3,4),已知数据见表2.9,试定位未知信号的位置坐标(x,y)
点击查看代码
from scipy.optimize import least_squares
import numpy as np
a=np.loadtxt('data2_47.txt')
x0=a[0]; y0=a[1]; d=a[2]
fx=lambda x: np.sqrt((x0-x[0])**2+(y0-x[1])**2)-d
s=least_squares(fx, np.random.rand(2))
print(s, '\n', '------------', '\n', s.x)
print("学号:3001")
例2.48求下列矩阵的最大模特征值和对应的特征向量
点击查看代码
from scipy.sparse.linalg import eigs
import numpy as np
a = np.array([[1, 2, 3], [2, 1, 3], [3, 3, 6]], dtype=float) #必须加float,否则出错
b, c = np.linalg.eig(a)
d, e = eigs(a, 1)
print('最大模特征值为:', d)
print('对应的特征向量为:\n', e)
print("学号:3001")
例2.49利用solve求下列符号代数方程的解ax**2+bx+c=0,其中x为未知数
点击查看代码
import sympy as sp
a, b, c, x=sp.symbols('a,b,c,x')
x0=sp.solve(a*x**2+b*x+c, x)
print(x0)
print("学号:3001")
例2.50求方程组的符号解
点击查看代码
#1
import sympy as sp
sp.var('x1,x2')
s=sp.solve([x1**2+x2**2-1,x1-x2],[x1,x2])
print(s)
print("学号:3001")
#2
import sympy as sp
x = sp.var('x:2') #定义符号数组
s = sp.solve([x[0]**2+x[1]**2-1,x[0]-x[1]], x)
print(s)
print("学号:3001")
例2.51求下列矩阵的特征值和特征向量的符号解
点击查看代码
import numpy as np
import sympy as sp
a = np.identity(4) #单位矩阵的另一种写法
b = np.rot90(a)
c = sp.Matrix(b)
print('特征值为:', c.eigenvals())
print('特征向量为:\n', c.eigenvects())
print("学号:3001")
例2.52已知店铺商品的销售量如表,画出商品销售趋势图
点击查看代码
import pandas as pd
import pylab as plt
plt.rc('font',family='SimHei') #用来正常显示中文标签
plt.rc('font',size=16) #设置显示字体大小
a=pd.read_excel("data2_52.xlsx", header=None)
b=a.values #提取其中的数据
x=b[0]; y=b[1:]
plt.plot(x,y[0],'-*b',label='铂金')
plt.plot(x,y[1],'--dr',label='铂金')
plt.xlabel('月份'); plt.ylabel('每月销量')
plt.legend(loc='upper left'); plt.grid(); plt.show()
print("学号:3001")
例2.53画出上表销售数据的柱状图
点击查看代码
import pandas as pd
import pylab as plt
plt.rc('font',family='SimHei') #用来正常显示中文标签
plt.rc('font',size=16) #设置显示字体大小
a=pd.read_excel("data2_52.xlsx",header=None)
b=a.T; b.plot(kind='bar'); plt.legend(['钻石', '铂金'])
plt.xticks(range(6), b[0], rotation=0); plt.show()
print("学号:3001")
例2.54把一个窗口分成3个子窗口,分别绘制如下3个子图:
一个柱状图 一个饼图 曲线y=sin(10x)/x
点击查看代码
import pylab as plt
import numpy as np
plt.style.use('default')
y1 = np.random.randint(2, 5, 6)
y1 = y1 / sum(y1)
plt.subplot(2, 2, 1)
str = ['Apple', 'grape', 'peach', 'pear', 'banana', 'pineapple']
plt.barh(str, y1)
plt.subplot(222)
plt.pie(y1, labels=str)
plt.subplot(212)
x2 = np.linspace(0.01, 10, 100)
y2 = np.sin(10*x2) / x2
plt.plot(x2, y2)
plt.xlabel('$x$')
plt.ylabel('$\\mathrm{sin}(10x)/x$')
plt.show()
print("学号:3001")
例2.55画出三维曲线x=s2sins,y=s2coss,z=s(s<-50,50>)的图形
点击查看代码
import pylab as plt
import numpy as np
ax=plt.axes(projection='3d') #设置三维图形模式
z=np.linspace(-50, 50, 1000)
x=z**2*np.sin(z); y=z**2*np.cos(z)
plt.plot(x, y, z, 'k'); plt.show()
print("学号:3001")
例2.56画出三维表面图z=50sin(x+y)
点击查看代码
import pylab as plt
import numpy as np
x=np.linspace(-4,4,100);
x,y=np.meshgrid(x,x)
z=50*np.sin(x+y);
ax=plt.axes(projection='3d')
ax.plot_surface(x, y, z, color='y')
plt.show()
print("学号:3001")
例2.57画出三维表面图z=sin(sqrt(x2+y2))
点击查看代码
import pylab as plt
import numpy as np
ax=plt.axes(projection='3d')
X = np.arange(-6, 6, 0.25)
Y = np.arange(-6, 6, 0.25)
X, Y = np.meshgrid(X, Y)
Z = np.sin(np.sqrt(X**2 + Y**2))
surf = ax.plot_surface(X, Y, Z, cmap='coolwarm')
plt.colorbar(surf); plt.show()
print("学号:3001")
第二章习题
习题2.1
点击查看代码
from cProfile import label
from re import X
import matplotlib.pyplot as plt
import numpy as np
from matplotlib import font_manager
my_font = font_manager.FontProperties(fname=
"C:/Windows/Fonts/msyh.ttc")
plt.figure(figsize=(6,5))#显示画布设置
plt.subplot(111)#第一个参数表示:行,第二个参数表示;列,第三个参数;当前图例中的激活位置
plt.xlabel(u'X数值',fontproperties=my_font)#x标签设置
plt.ylabel(u'Y数值',fontproperties=my_font)#y标签设置
plt.title(u"函数图像",fontproperties=my_font,fontsize=16)#图像名称设置
x1=np.linspace(-10,10,100)#-10,10均分100份
e= 2.718281828459
#绘制chx双曲余弦
y1=((e**x1)+(e**-x1))/2
plt.plot(x1,y1,label=u'chx')
#绘制双曲正弦shx
y2=((e**x1)-(e**-x1))/2
plt.plot(x1,y2,label=u'shx')
#绘制,1/2e^x
y3=1/2*(e**x1)
plt.plot(x1,y3,label=u'1/2e^x')
plt.legend(prop=my_font)#显示函数文本标注。
plt.show()
print("学号:3001")
习题2.2
点击查看代码
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad
def fun(t, x):
return np.exp(-t) * (t ** (x - 1))
x = np.linspace(0, 10, 100) # x 的范围
y = [quad(fun, 0, np.inf, args=i)[0] for i in x] # 计算积分
plt.plot(x, y)
plt.xlabel('x')
plt.ylabel('$ y = \int_0^{\infty} e^{-t} \cdot t^{x-1} dt $')
plt.grid(True)
plt.show()
print("学号:3001")
习题2.3
点击查看代码
import numpy as np
import matplotlib.pyplot as plt
k_values = [1, 2, 3, 4, 5, 6] # k的取值
x = np.linspace(-10, 10, 100) # x的范围
for k in k_values:
y = k * x ** 2 + 2 * k # 计算y值
label = f'k={k}' # 设置标注
plt.plot(x, y, label=label) # 绘制曲线
plt.xlabel('x') # 添加x轴标签
plt.ylabel('y') # 添加y轴标签
plt.legend() # 添加图例
plt.grid(True) # 添加网格线
plt.show() # 显示图形
print("学号:3001")
习题2.4
点击查看代码
import numpy as np
import matplotlib.pyplot as plt
plt.rc('font', family='SimHei') # 用来正常显示中文标签
plt.rc('axes', unicode_minus=False) # 用来正常显示负号
k_values = [1, 2, 3, 4, 5, 6] # k的取值
x = np.linspace(-10, 10, 100) # x的范围
fig, axs = plt.subplots(2, 3, figsize=(10, 6)) # 创建2行3列的子窗口
for i, k in enumerate(k_values):
y = k * x ** 2 + 2 * k # 计算y值
row = i // 3 # 计算子窗口所在的行数
col = i % 3 # 计算子窗口所在的列数
ax = axs[row, col] # 获取当前子窗口
label = f'k={k}' # 设置标注
ax.plot(x, y, label=label) # 绘制曲线
ax.set_xlabel('x') # 添加x轴标签
ax.set_ylabel('y') # 添加y轴标签
ax.set_title(f'图 {i+1}') # 添加子窗口标题
ax.legend() # 添加图例
ax.grid(True) # 添加网格线
plt.tight_layout() # 自动调整子窗口布局
plt.show() # 显示图形
print("学号:3001")
习题2.5
(1)
点击查看代码
import numpy as np
import matplotlib.pyplot as plt
a = 2
b = np.sqrt(10)
c = np.sqrt(8)
phi = np.arange(0, 2*np.pi+0.1, 0.1)
theta = np.arange(-1, 1.1, 0.1)[:, np.newaxis]
x = a * np.cosh(theta) * np.cos(phi)
y = b * np.cosh(theta) * np.sin(phi)
z = c * np.sinh(theta)
fig = plt.figure()
ax = fig.add_subplot(projection='3d')
ax.plot_surface(x, y, z, cmap='viridis')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
ax.set_box_aspect([1, 1, 1])
ax.set_title('$\\frac{x^2}{4}+\\frac{y^2}{10}-\\frac{z^2}{8}=1$')
plt.show()
print("学号:3001")
(2)
点击查看代码
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# 创建一个三维坐标系
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
# 生成x,y的网格
x = np.linspace(-50, 50, 1000)
y = np.linspace(-50, 50, 1000)
X, Y = np.meshgrid(x, y)
# 计算z的值
Z = (X**2)/4+(Y**2)/6
# 绘制二次曲面
ax.plot_surface(X, Y, Z, cmap='viridis')
# 添加坐标轴标签
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
# 添加标题
ax.set_title('$(\\frac{x^2}{4}+\\frac{y^2}{6})=z$')
# 设置z轴的范围
ax.set_zlim(0, 1000)
# 显示图形
plt.show()
print("学号:3001")
习题2.6
点击查看代码
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# 模拟高程数据(假设数据已经过某种方式插值或生成)
# 这里我们创建一个简单的40x50网格,并填充随机高程值
x = np.linspace(0, 43.65, 40)
y = np.linspace(0, 58.2, 50)
X, Y = np.meshgrid(x, y)
Z = np.sin(np.sqrt(X**2 + Y**2)) * 100 # 使用一个简单的函数来生成高程数据
# 绘制三维表面图
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(121, projection='3d')
ax.plot_surface(X, Y, Z, cmap='viridis', edgecolor='none')
ax.set_xlabel('X (km)')
ax.set_ylabel('Y (km)')
ax.set_zlabel('Elevation (m)')
ax.set_title('3D Surface Plot of Elevation Data')
# 绘制等高线图
plt.subplot(122)
CS = plt.contour(X, Y, Z, colors='k')
plt.clabel(CS, inline=1, fontsize=10)
# 标注点 A(30,0) 和 B(43,30)
# 注意:由于X和Y是网格坐标,我们需要找到最接近这些值的索引
idx_a_x = np.argmin(np.abs(x - 30))
idx_a_y = np.argmin(np.abs(y - 0))
idx_b_x = np.argmin(np.abs(x - 43))
idx_b_y = np.argmin(np.abs(y - 30))
plt.plot(x[idx_a_x], y[idx_a_y], 'ro', markersize=5, label='A(30,0)')
plt.plot(x[idx_b_x], y[idx_b_y], 'go', markersize=5, label='B(43,30)')
plt.xlabel('X (km)')
plt.ylabel('Y (km)')
plt.title('Contour Plot of Elevation Data with Points A and B')
plt.legend()
# 计算地表面积的近似值(忽略地形起伏)
real_area = 43.65 * 58.2
print(f"Actual Surface Area (ignoring elevation changes): {real_area} km^2")
# 显示图形
plt.tight_layout()
plt.show()
print("学号:3001")
习题2.7
点击查看代码
import numpy as np
# 定义系数矩阵A和常数项向量b
A = np.array([[4, 2, -1],
[3, -1, 2],
[11, 3, 0]])
b = np.array([2, 10, 8])
# 使用numpy的lstsq求解最小二乘解
x, residuals, rank, s = np.linalg.lstsq(A, b, rcond=None)
print("最小二乘解为:")
print(x)
# 打印残差和矩阵A的秩
print("残差为:", residuals)
print("矩阵A的秩为:", rank)
print("学号:3001")
print("\n")
import numpy as np
# 定义系数矩阵A和常数项向量b
A = np.array([[2, 3, 1],
[1, -2, 4],
[3, 8, -2],
[4, -1, 9]])
b = np.array([4, -5, 13, -6])
# 使用numpy的lstsq函数求解最小二乘解
# 对于这个特定的问题,由于方程数和未知数数量相同,且没有矛盾,lstsq将给出唯一解
x, residuals, rank, s = np.linalg.lstsq(A, b, rcond=None)
# 输出解
print("解 x:", x)
# 验证解是否正确(可选)
print("验证结果:", np.dot(A, x))
# 检查是否精确等于b(对于精确解,这应该非常接近)
print("与b的误差:", np.linalg.norm(np.dot(A, x) - b))
# 计算系数矩阵的秩(可选,以确认方程组是否有唯一解)
print("系数矩阵的秩:", np.linalg.matrix_rank(A))
# 由于秩等于未知数数量,且没有矛盾,我们可以确信有一个唯一解
print("学号:3001")
习题2.8
点击查看代码
import numpy as np
# 生成系数矩阵A
A = np.zeros((1000, 1000))
np.fill_diagonal(A, 4)
np.fill_diagonal(A[:, 1:], 1)
np.fill_diagonal(A[1:, :], 1)
# 生成常数向量b
b = np.arange(1, 1001)
# 判断解的情况
if np.linalg.matrix_rank(A) == np.linalg.matrix_rank(np.column_stack((A, b))):
if np.linalg.matrix_rank(A) == A.shape[1]:
print("线性方程组有唯一解")
x_unique = np.linalg.solve(A, b)
print("唯一解 x =", x_unique)
else:
print("线性方程组有无穷多解")
x_least_squares = np.linalg.lstsq(A, b, rcond=None)[0]
print("最小二乘解 x =", x_least_squares)
else:
print("线性方程组无解")
x_least_squares = np.linalg.lstsq(A, b, rcond=None)[0]
print("最小二乘解 x =", x_least_squares)
x_min_norm = np.linalg.pinv(A).dot(b)
print("最小范数解 x =", x_min_norm)
import numpy as np
# 生成系数矩阵A
A = np.zeros((1000, 1000))
np.fill_diagonal(A, 4)
np.fill_diagonal(A[:, 1:], 1)
np.fill_diagonal(A[1:, :], 1)
# 生成常数向量b
b = np.arange(1, 1001)
# 判断解的情况
if np.linalg.matrix_rank(A) == np.linalg.matrix_rank(np.column_stack((A, b))):
if np.linalg.matrix_rank(A) == A.shape[1]:
print("线性方程组有唯一解")
x_unique = np.linalg.solve(A, b)
print("唯一解 x =", x_unique)
else:
print("线性方程组有无穷多解")
x_least_squares = np.linalg.lstsq(A, b, rcond=None)[0]
print("最小二乘解 x =", x_least_squares)
else:
print("线性方程组无解")
x_least_squares = np.linalg.lstsq(A, b, rcond=None)[0]
print("最小二乘解 x =", x_least_squares)
x_min_norm = np.linalg.pinv(A).dot(b)
print("最小范数解 x =", x_min_norm)
print("学号:3001")
习题2.9
点击查看代码
from sympy import symbols, Eq, solve
# 定义未知数
x, y = symbols('x y')
# 定义非线性方程组
equations = (Eq(x**2 - y - x - 3, 0), Eq(x + 3*y - 2, 0))
# 求解符号解
symbolic_solution = solve(equations, (x, y))
print("符号解:", symbolic_solution)
# 数值解需要使用数值方法求解
from scipy.optimize import fsolve
# 定义非线性方程组的函数
def equations(variables):
x, y = variables
return [x**2 - y - x - 3, x + 3*y - 2]
# 初始猜测值
initial_guess = [1, 1]
# 求解数值解
numeric_solution = fsolve(equations, initial_guess)
print("数值解:", numeric_solution)
print("学号:3001")
习题2.10
点击查看代码
from scipy.integrate import quad
import numpy as np
#抛物线旋转体
def V1_quad(y):
return np.pi * (4*y - y**2)
V1_corrected, _ = quad(V1_quad, 1, 3)
V2 = 0.5 * (4/3) * np.pi * 2**3 - (1/3) * np.pi * 2**2 * 1
# 计算总体积
total_volume_corrected = V1_corrected + V2
print("体积:",total_volume_corrected)
import numpy as np
import math
# 圆柱面部分
V2 = 4 * math.pi # 体积
y2 = 0.5 # 重心y坐标
# 假设水的密度 rho = 1000 kg/m^3
rho = 1000 # 水的密度
g = 9.81 # 重力加速度
# 计算重力势能变化和所需功
# 假设水被抽到无穷高(实际中可能有限制),这里以 y = 10(远大于容器高度)为例
final_y = 10
delta_E_p = rho * V2 * g * (final_y - y2)
W = delta_E_p
print(f"圆柱面部分所需功: {W} 焦耳")
print("学号:3001")
习题2.11
点击查看代码
import numpy as np
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)
def equations(variables):
x1, x2, y1, y2 = variables[:4]
eq1 = 2 * x1 - (3 * f(y1) + 4 * g(y2) - 1)
eq2 = 3 * x2 - (2 * f(y1) + 6 * g(y2) - 2)
eq3 = y1 - (f(x1) + 3 * g(x2) - 3)
eq4 = 5 * y2 - (4 * f(x1) + 6 * g(x2) - 1)
return [eq1, eq2, eq3, eq4]
initial_guess = [0, 0, 0, 0] # 初始猜测值
numeric_solution = fsolve(equations, initial_guess)
print("数值解:", numeric_solution)
print("学号:3001")
习题2.12
点击查看代码
import numpy as np
from scipy.linalg import eig
# 定义矩阵
A = np.array([[-1, 1, 0],
[-4, 3, 0],
[1, 0, 2]])
# 计算特征值和特征向量
eigenvalues, eigenvectors = eig(A)
# 打印特征值
print("特征值:")
print(eigenvalues)
# 打印特征向量
print("特征向量:")
for i in range(eigenvectors.shape[1]):
print(f"特征值 {eigenvalues[i]:.2f} 对应的特征向量:")
print(eigenvectors[:, i].real) # 取实部,因为有时特征向量会有复数部分,但在这个特定例子中它们是实数
print("学号:3001")
习题2.13
点击查看代码
import numpy as np
from scipy.optimize import least_squares
def f(x):
return (np.abs(x + 1) - np.abs(x - 1)) / 2 + np.sin(x)
def g(x):
return (np.abs(x + 3) - np.abs(x - 3)) / 2 + np.cos(x)
def equations(variables):
x1, x2, y1, y2 = variables[:4]
eq1 = 2 * x1 - (3 * f(y1) + 4 * g(y2) - 1)
eq2 = 3 * x2 - (2 * f(y1) + 6 * g(y2) - 2)
eq3 = y1 - (f(x1) + 3 * g(x2) - 3)
eq4 = 5 * y2 - (4 * f(x1) + 6 * g(x2) - 1)
eq5 = x1 + y1 - (f(y2) + g(x2) - 2)
eq6 = x2 - 3 * y2 - (2 * f(x1) - 10 * g(y1) - 5)
return [eq1, eq2, eq3, eq4, eq5, eq6]
initial_guess = [0, 0, 0, 0] # 初始猜测值
result = least_squares(equations, initial_guess)
numeric_solution = result.x
print("最小二乘解:", numeric_solution)
print("学号:3001")