python数学建模第八章

8.4

from scipy.integrate import odeint
import numpy as np
import pylab as plt

plt.rc('font',size=15); plt.rc('text',usetex=True)
dz=lambda z,t: [-z[0]**3-z[1],z[0]-z[1]**3]
t0=np.linspace(0,30,201)
s=odeint(dz,[1,0.5],t0)
plt.subplot(121)
plt.plot(t0,s[:,0],'--',label='$x(t)$')
plt.plot(t0,s[:,1],label='$y(t)$')
plt.legend(); plt.subplot(122)
plt.plot(s[:,0],s[:,1]); plt.show()

8.5

from scipy.integrate import odeint
import numpy as np
import pylab as plt

plt.rc('font',size=15); plt.rc('text',usetex=True)
dy=lambda y,t: [y[1],y[2],-3*y[0]*y[2]+2*y[1]**2-y[3],y[4],-2.1*y[0]*y[4]]
t0=np.linspace(0,10,51); y0=[0,0,0.68,1,-0.5]
s=odeint(dy,y0,t0)
plt.plot(t0,s[:,0],'--*',label='$f(\eta)$')
plt.plot(t0,s[:,3],'p-',label='$T(\eta)$')
plt.xlabel('$\eta$'); plt.legend(); plt.show()

8.8

import numpy as np
import pylab as plt

a = 1 - 0.2**(1/12); m = 1.109 * 10 ** 5
w3 = 17.86; w4 = 22.99
X=[]; Z=[]; N=[]
for k in np.arange(0, 0.875, 0.001):
    x1 = 1.22*10**11*(1-1/(m*(1-a-0.42*k)**8*(1-a)**24*
        (1/2+(1-a-k)**8*(1-a)**4/(1-(1-a-k)**8*(1-a)**4))))
    x2 = (1-a)**12*x1; x3 = (1-a)**12*x2;
    x4 = (1-a-0.42*k)**8*(1-a)**4/(1-(1-a-k)**8*(1-a)**4)*x3
    X.append([x1,x2,x3,x4])
    n = (m*(1-a-0.42*k)**8*(1-a)**24*(1/2+(1-a-k)**8*(1-a)**4/
        (1-(1-a-k)**8*(1-a)**4))*x1)
    N.append(n)
    z = (0.42*k*w3*(1-(1-a-0.42*k)**8)/(a+0.42*k)*x3+
         k*w4*(1-(1-a-k)**8)/(a+k)*x4)
    Z.append(z)
mz = max(np.array(Z)); ind = np.argmax(Z)
k4 = 0.001*ind; k3 = 0.42 * k4  #最优捕捞强度
print('最大生产量：', round(mz,4))
print('各年龄组鱼群数：', np.round(X[ind],4))
print('k4=',round(k4,4)); print('k3=',round(k3,4))
plt.rc('text', usetex=True); plt.rc('font', size=16)
plt.plot(np.arange(0, 0.875, 0.001), Z)
plt.ylabel('$z$', rotation=0)
plt.xlabel('$k$'); plt.show()

8.9

r = 0.0036; N = 360; Q = 400000
x = round((1+r)**N*Q*r/((1+r)**N-1),2)
xt = x*N  #总还款额
print('x=',x); print('xt=',xt)