파이썬 질문

파이썬으로 이송확산 방정식을 통한 수치해석에 대해 공부하고 있는 학생입니다.. 자꾸 코드를 실행하면 c와 rhs배열에 nan이 들어가는데 그 이유를 모르겠습니다..

#Library import import numpy as np import math import matplotlib.pyplot as plt

#Definitions and initializations of variable arrays m = 100 v = 1.0 dx = 1.0 dt = 0.02 nmax = 2250 # T / dt (T=45) x = np.zeros(shape = m+1) # 0~100 array  c = np.zeros(shape = m+1) rhs = np.zeros(shape = m+1) # => 행벡터 a1 = np.zeros(shape = m+1) b1 = np.zeros(shape = m+1) c1 = np.zeros(shape = m+1) d = 0.0 e = 0.0 f = 0.0 c0 = np.zeros(shape = m+1) ce = np.zeros(shape = m+1)

# Selection of scheme alpha = float(input('input alpha : ')) #키보드 입력 , cf) 강제 형변환  theta = float (input('input theta : ')) #0, 0.5, 1.0  (theta가 1일때 implicit method, 0일때 explicit method) pe = float(input('input pe : '))#0.1 , 1 scheme = float(input('input scheme :')) 

#Grid generation #% when i = 0, its distance is at (1*dx)m from origin for i in range (0, m+1):     x[i] = i*dx      #initial condition for i in range (0, m+1):     if i == 0:         c[i] = 1.0     else :         c[i] = 0.0

#Save the initial condition for i in range (0, m+1):     c0[i]=c[i]         #========================================== #Main Program #========================================== #Implicit method  #Time integral

for iter in range (1,nmax+1): #for 1문     if iter == round(iter/100)*100: #100번 iteration마다 찍고 싶을 때 사용         print('{0}'.format(iter))     for i in range (1,m): ## for 2문         a1[i] = (-v*alpha*theta/(2*dx))-v*(1-alpha)*theta/dx-v*theta/(pe*dx);         b1[i] = 1/dt + (v*theta*(1-alpha)/dx) + (2*v*theta/(pe*dx))         c1[i] = (v*alpha*theta/(2*dx))-(v*theta/(pe*dx))         #RHS Terms          d = (v*alpha*(1-theta)/(2*dx))+(v*(1-alpha)*(1-theta)/dx)+(v*(1-theta)/(pe*dx))         e = 1/dt - (v*(1-alpha)*(1-theta)/dx)-(2*v*(1-theta)/(pe*dx))         f = (-v*alpha*(1-theta)/(2*dx)) + (v*(1-theta)/(pe*dx))         rhs[i] = d*c[i-1]+e*c[i]+f*c[i+1] ##end of for 2문              rhs[1]=rhs[1]-a1[1]*c[0]     rhs[m-1]=rhs[m-1]-c1[m-1]*c[m]           #Thomas algorithm     for i in range(2,m): ##for 3문         r=a1[i]/b1[i-1]         b1[i]=b1[i]-r*c1[i-1]         rhs[i]=rhs[i]-r*rhs[i-1] ##end of for 3문           rhs[m-1]=rhs[m-1]/b1[m-1];     for i in range(m-2,0,-1): ##for 4문         rhs[i]=(rhs[i]-c1[i]*rhs[i+1]/b1[i])##end of for 4문     #end of Thomas algorithm          #renewal of concentration at the next time step     for i in range (1, m): ##for 5문         c[i]=rhs[i] ## end of for 5문     #renwal of boundary conditions      c[0] = 1.0 #continous source injection     c[m] = 0.0 #Neumann boundary condition     # end of for 1문      #Exact Solution for i in range (0, m+1):     d = v*dx/pe     y1 = (x[i]-v*dt*nmax)/(2*math.sqrt(d*dt*nmax))     y2 = (x[i]+v*dt*nmax)/(2*math.sqrt(d*dt*nmax))     ce[i]=c[0]/2*(math.erfc(y1)+math.exp(v*x[i]/d)*math.erfc(y2))