python数学实验与建模（习题3）

3.1

# In[3.1]
from sympy import integrate,sqrt,exp,oo,sin
from sympy.abc import x

print(integrate(sqrt(1+4*x),(x,0,1)))
print(integrate(exp(-x)*sin(x),(x,0,oo)))

3.2

# In[3.2]
from sympy import *
from sympy.abc import x
# 符号解
s=solve(x**3-4*x**2+6*x-8,x)
print(s)
# 数值解
print(s[0].n())

3.3

# In[3.3]
from sympy import *
from sympy.abc import x,y
# 符号解
s=solve([x**2-y-x-3,x+3*y-6],[x,y])
print(solve([x**2-y-x-3,x+3*y-6],[x,y]))
# 数值解
print(s[i].n() for i in range(len(s)))

3.4

# In[3.4]
from sympy import *
from sympy.abc import x,y

y=Function('y')

eq=diff(y(x),x,2)+diff(y(x),x,1)-x*cos(2*x)

print(dsolve(eq,y(x),ics={y(0):1,y(2):3}))

3.5

# In[3.5]
import sympy as sp
import numpy as np
A1=np.arange(1,7).reshape(3,2)
A2=np.array([[1,1],
            [2,2],
            [3,4]])

A3=np.array([[2,6]])
A4=np.array([[3,2]])

M1 = np.concatenate((A1, A2), axis=1)
M2 = np.concatenate((A3, A4), axis=1)
M = np.concatenate((M1, M2), axis=0)

print(M)
print(sp.det(sp.Matrix(M)))

3.6

3.6.1

# In[3.6.1]
# 齐次线性方程组 求基础解系
import sympy as sp
import numpy as np
A=sp.Matrix([[1,2,1,-1],[3,6,-1,-3],[5,10,1,-5]])
b=sp.Matrix([0,0,0,0])
b.transpose()
print(A.nullspace())

3.6.2

# In[3.6.2]
# 非齐次线性方程组 求通解
import sympy as sp
import numpy as np
A=sp.Matrix([[2,1,-1,1],[4,2,-2,1],[2,1,-1,-1]])
b=sp.Matrix([1,2,1])
b.transpose()
C=np.hstack([A,b])
C=sp.Matrix(C)
D=C[0:2,0:5]
print("方程通解为：\n",D.nullspace())

3.7

3.7.1

# In[3.7.1]
import sympy as sp
import numpy as np
from scipy.optimize import fsolve
A=sp.Matrix([[4,2,-1],[3,-1,2],[11,3,0]])
b=sp.Matrix([2,10,8])
b.transpose()
C=np.hstack([A,b])
D=sp.Matrix(C)
# 确定问题为求解非线性方程组
print(D.rref())

f=lambda x:[4*x[0]+2*x[2]-x[2]-2,3*x[0]-x[1]+2*x[2]-10,11*x[0]+3*x[1]-8]
print(fsolve(f,[1,1,1]))

3.7.2

# In[3.7.2]
import sympy as sp
import numpy as np
from scipy.optimize import fsolve
A=sp.Matrix([[2,3,1],[1,-2,4],[3,8,-2],[4,-1,9]])
b=sp.Matrix([4,-5,13,-6])
b.transpose()
C=np.hstack([A,b])
D=sp.Matrix(C)
# 确定问题为求解非齐次线性方程组
print(D.rref()) 
E=D[0:2,0:5]
print("方程通解为：\n",D.nullspace())

3.8

# In[3.8]
import sympy as sp
A=sp.Matrix([[6,2,4],[2,3,2],[4,2,6]])
print("A的特征值为：",A.eigenvals())
print("A的特征向量为：",A.eigenvects())

3.10

from sympy import *
from sympy.abc import x,y
from sympy.plotting import plot
# 定义函数
y=cos((x**2+1)**0.5)
# 间隔取不同展开阶数
for k in range(1,6,2):
    # 写展开式
    print(y.series(x,0,k))
    # 画展开图
plot(y,series(y,x,0,1).removeO(),series(y,x,0,3).removeO(),series(y,x,0,5).removeO())

 

3.11

# In[3.11]
import numpy as np
import matplotlib.pyplot as plt
N=2;vr=3;vd=4.5
xy=np.array([[0,0],[20,0]])
xyn=np.zeros((2,2))
T=np.linspace(0,7,200)
dt=T[1]-T[0]
Txy=xy
for n in range(0,len(T)):
        xy[1]=[20,vr*dt*n]
        dxy=xy[1]-xy[0]
        dd=dxy/np.linalg.norm(dxy)
        xyn[0]=xy[0]+dd*vd*dt
        xyn[1]=xy[1]
        Txy=np.c_[Txy,xyn]
        xy=xyn.copy()
for i in [0,1]:
    plt.plot(Txy[i,::2],Txy[i,1::2])

可以改变代码中T=np.linspace(0,7,200)中第二个参数，并不断绘制图像，上图是时间终点取7个单位时间时的结果，由此估计需要七秒左右狐狸追上兔子

3.12

 # In[3.12]
from sympy import integrate,sqrt,exp,oo,sin
from sympy.abc import x
from scipy.integrate import dblquad,tplquad
import numpy as np
# 1
print("I1:",(integrate((exp(-x)*sin((x**2+2)**0.5)),(x,0,oo))).n())
# 2
# 注意先积分y后积分x
f1=lambda y,x:x**2+2*y**2
bd1=lambda x:np.sqrt(x)
bd2=lambda x:x-2

f=(dblquad(f1,1,4,lambda x:bd2(x),lambda x:bd1(x))+
              dblquad(f1,0,1,lambda x:-bd1(x),lambda x:bd1(x)))
print("I2:",f[0]+f[2])
# 3
# 体会两种不同的写法 在设置匿名函数的时候需要注意传参个数
f2=lambda theta,r,z:r*z
t=lambda z:np.sqrt(z)
print('I3:',tplquad(f2,0,4,0,lambda z:t(z),0,2*np.pi))

f2=lambda r,theta,z:r*z
t=lambda z,theta:np.sqrt(z)
print('I3:',tplquad(f2,0,4,0,2*np.pi,0,lambda z,theta:t(z,theta)))

3.13

# In[3.13]
from numpy import exp,sin
from scipy.optimize import fminbound
# lambda函数要求输入数字而非符号 要充分理解符号解和数值解的含义
y=lambda x:2*exp(-x)*sin(x)
t=lambda x:-2*exp(-x)*sin(x)
x1=fminbound(y,0,3)
x2=fminbound(t,0,3)
print("极小点：{},极大点：{}".format(x1,x2))