text: import math
  import numpy as np
  import matplotlib.pyplot as plt
  
  
  def func(x,k):
      #k = 0.707 阻尼比系数
      b = math.acos(k) #阻尼角
      wn = 1 #无阻尼振荡角频率
      c1 = math.sqrt(1-k**2)  
      wd = wn*c1 #阻尼振荡角频率
      ans = 1-math.exp(-k*wn*x)*math.sin(wd*x+b)/c1
      return ans
  
  def realTime(k):
      t = np.linspace(0,30,num=30000) # delta t = 0.001
      seqk = [k]*30000
      y = list(map(func,t,seqk))
      realtime = 0
      ttt = []
      for i in range(1,30000+1):
          if  y[-i]<= 1-0.05 or y[-i] >= 1+0.05:
              realtime=t[-i]
              break
      return realtime
  
  def envelopeTime(k):
      return -(math.log(0.05)+math.log(math.sqrt(1-k**2)))/k
  def estimateTime(k):
      return 3/k
  
  def solve():
      k = np.linspace(0.1,0.999,num=1000)
      realtime = list(map(realTime,k))
      envelopetime = list(map(envelopeTime,k))
      estimatetime = list(map(estimateTime,k))
      plt.figure()
      plt.plot(k,realtime,label='Real Time')
      plt.plot(k,envelopetime,label ='Envelope Time')
      plt.plot(k,estimatetime,label = 'Estimate Time')
      plt.legend()
      plt.xlabel(r'$\zeta$')
      plt.ylabel(r'$\omega_nt$')
      plt.show()
      
  
  solve()
  