I have the following program transcribed in Python from MATLAB
import numpy as np
from math import exp,ceil
def initial (x):
init=0.09*np.exp(-x**2/50)
return init
a = -20
b= 20
T=4
M=500
rhom=0.2
n=10
drho=rhom/n
rho=np.linspace(0,rhom,n+1)#lim inf, lim sup, "step"
def flux(rho):
vf=15
rhomax=0.2
fl=vf*(1-rho/rhomax)*rho
return fl
frho=flux(rho)
qmax=max(frho)#Busca el máximo de la lista
s=np.argmax(frho)#devuelve el índice del valor máximo
rhostar=rho.item(s)#devuelce el elemento en la posicion "i"
#CFL Condition
lamda=max(abs(frho[1:n]-frho[0:n-1]))/drho
h=(b-a)/M
ka=0.5*h/lamda
N=ceil(T/ka)
k=T/N
#Inicializar:
kio=np.insert(np.arange(a+h/2,b-h/2,h),len(np.arange(a+h/2,b-h/2,h)),b-h/2)# np.insert(lista,length,valor_a_agragar)
xticks=kio.transpose()
U=np.zeros((M,N+1))
U[:,0]=initial(xticks)
for j in range(N):
rhol=U[0:M-1,j-1]#U(1:M-1,j)
rhor=U[1:M+1,j-1]#U(2:M,j)
qval=flux(U[:,j-1])
ql=qval[0:M-1]
qr=qval[1:M+1]
case1=(rhol<=rhor)
case2=(rhol>rhostar)&(rhostar>rhor)
case3= ~(case1|case2)
Q1=np.sum(case1)*np.fmin(ql,qr)+np.sum(case2)*qmax+np.sum(case3)*np.fmax(ql,qr)
Q2=np.insert(Q1,len(Q1),qval[M-1])
Q=np.insert(Q2,0,flux(initial(a)))#np.insert("array","indice donde quiero que vaya el elemento","elemento")
U[:,j+1]=U[:,j]+(k/h)*(Q[0:M]-Q[1:M+1])
print(U)
When I run the program I get the following
Warning (from warnings module):
File "C:\Users\ALEX EDUARDO\Documents\HolaMundoPython.py", line 20
fl=vf*(1-rho/rhomax)*rho
RuntimeWarning: overflow encountered in multiply
Warning (from warnings module):
File "C:\Users\ALEX EDUARDO\Documents\HolaMundoPython.py", line 50
U[:,j+1]=U[:,j]+(k/h)*(Q[0:M]-Q[1:M+1])
RuntimeWarning: invalid value encountered in subtract
Warning (from warnings module):
File "C:\Users\ALEX EDUARDO\Documents\HolaMundoPython.py", line 47
Q1=np.sum(case1)*np.fmin(ql,qr)+np.sum(case2)*qmax+np.sum(case3)*np.fmax(ql,qr)
RuntimeWarning: overflow encountered in multiply
Warning (from warnings module):
File "C:\Users\ALEX EDUARDO\Documents\HolaMundoPython.py", line 44
case1=(rhol<=rhor)
RuntimeWarning: invalid value encountered in less_equal
Warning (from warnings module):
File "C:\Users\ALEX EDUARDO\Documents\HolaMundoPython.py", line 45
case2=(rhol>rhostar)&(rhostar>rhor)
RuntimeWarning: invalid value encountered in greater
Warning (from warnings module):
File "C:\Users\ALEX EDUARDO\Documents\HolaMundoPython.py", line 45
case2=(rhol>rhostar)&(rhostar>rhor)
RuntimeWarning: invalid value encountered in less
Warning (from warnings module):
File "C:\Users\ALEX EDUARDO\Documents\HolaMundoPython.py", line 47
Q1=np.sum(case1)*np.fmin(ql,qr)+np.sum(case2)*qmax+np.sum(case3)*np.fmax(ql,qr)
RuntimeWarning: invalid value encountered in multiply
[[ 3.11723957e-05 4.79429950e-05 -8.85935460e-03 ... nan
nan nan]
[ 3.32241476e-05 3.32241476e-05 -5.52799370e-04 ... nan
nan nan]
[ 3.54018808e-05 3.54018808e-05 -5.86503957e-04 ... nan
nan nan]
...
[ 3.54018808e-05 3.54018808e-05 6.57307718e-04 ... nan
nan nan]
[ 3.32241476e-05 3.32241476e-05 6.19247665e-04 ... nan
nan nan]
[ 3.11723957e-05 3.11723957e-05 8.93792529e-03 ... nan
nan nan]]
And the result of the matrix U should be as in the image (MATLAB):
