preheting

chi4_v2.py

from numpy import zeros, linspace, pi, cos,sin, array
import numpy as np
import math
import matplotlib.pyplot as plt

def getCalc(order, w, Phi, phiOrder,phim, T, dt):
  print("order="+str(order)+" w="+str(w)+" Phi="+str(Phi))

  # omega = 2
  # m_phi = 50
  # Period = 2*pi/m_phi
#   dt = T / 1000
  # T = 100*Period
  N_t = int(round(T/dt))
  t = linspace(0, N_t*dt, N_t + 1)
  # Lambda = 2
  # Phi = 10



  u = zeros(N_t+1)
  v = zeros(N_t+1)
  lognum = zeros(N_t+1)

  # 初期条件
  X_0 = 1
  V_0 = 0
  u[0] = X_0
  v[0] = V_0
  lognum[0] = 0
  maxlognum = 0
  maxU = 1

  # 差分計算
  for n in range(N_t):
    u[n+1] = u[n] + dt*v[n]
    v[n+1] = v[n] - dt*w*u[n+1] - dt*(Phi*(sin(phim*(n+1)*dt)))**phiOrder*u[n+1]**order
    # v[n+1] = v[n] - dt*w*u[n+1] - dt*(Phi*5)**phiOrder*u[n+1]**order

    lognum[n+1] = math.log(((math.sqrt(w))/2)*((v[n+1]**2)/(w**2) + u[n+1]**2))
    if maxlognum < lognum[n+1]:
      maxlognum = lognum[n+1]
    if maxU < u[n+1]:
      maxU = u[n+1]
  
  print("last log(n)="+str(lognum[N_t-1]))
  print("max log(n)="+str(maxlognum))

  return t,u,v,lognum,maxU,maxlognum

def use_fowardeuler2():  

  t = []
  ulist = []
  u_rescale_list = []
  vlist = []
  lognumlist = []
  maxUlist = []
  maxlognumlist = []

  # y'' - (a - q cos(2t)) y^order == 0
  maxnum_m = 1 # max order
  maxnum = 1 #  parameter
#   T = 28 # maxtime
  T = 23
  dt = 0.0001
#  q = [0.01*n for n in range(maxnum)] # narrow resonance
  Phi = [200 for n in range(maxnum)] # broad resonance
  phiOrder = [2 for n in range(maxnum)]
  w = [1 for n in range(maxnum)]
  phim = [ math.pi  for n in range(maxnum)]
#   a = [100 for n in range(maxnum)]
  order = [3 for n in range(maxnum_m)]
#   alpha = [-0.000000000000001 for n in range(maxnum)]

#   order = [1+2*n for n in range(maxnum_m)]

  for m in range(maxnum_m):
    for n in range(maxnum):
      t,u,v,lognum,maxU,maxlognum = getCalc(order[m],w[n],Phi[n],phiOrder[n],phim[n], T, dt)
      u_rescale = u * maxlognum/maxU
      ulist.append(u)
      vlist.append(v)
      u_rescale_list.append(u_rescale)
      lognumlist.append(lognum)
      maxUlist.append(maxU)
      maxlognumlist.append(maxlognum)

  fig = plt.figure()


  for m in range(maxnum_m):
    for n in range(maxnum):
    #   plt.plot(t, ulist[m*maxnum+n], label="n="+str(order[m])+" "+"a="+str(a[n])+" q="+str(q[n]))   
      plt.plot(t, lognumlist[m*maxnum+n],label="w="+str(w[n])+" "+"Phi="+str(Phi[n])+" n="+str(order[m]))    

    #  fig.legend(l, str(n), "upper left")
  # l1,l2 = plt.plot(t,u_rescale, "b-", t, lognum, "r--")
  # l1 = plt.plot(t,u, "b-")
  # fig.legend((l1), ("numerical"), "upper")

  # fig.legend((l1,l2), ("numerical", "lognum"), "upper left")
  fig.legend()
  plt.xlabel("t")
  plt.ylabel("log(n)")
  plt.show()

use_fowardeuler2()