RCIP

CIPは非常に優秀な移流方程式に対する解法ですが、不連続部に若干のオーバーシュート、アンダーシュートが残ります。

これは補間関数に3次関数を選択したためであり、これを避けるために、別の関数を使えば良いじゃん！というのが、RCIP(Rational CIP)基本的な考え方です。有理関数を使えばオーバーシュートアンダーシュートを抑えられます。

ソースを晒しておきます。まあ、汚いです。

#advection #################################################

# RCIP ########### ########################################

############################################################

#import#####################################################

from pylab import * 

#set parameters ############################################

in_ff = open('test.txt')

imx=100

dx=1.0

dt=0.2

u=0.2

alpha = 1.0

eps = 1.0e-15

print "-"*40

print "u  =",u

print "dx =",dx

print "dt =",dt

print "-"*40

############################################################

#initial condition ################################

f= [0]

fn=  [0]

df = [0]

dfn = [0]

fini=[0]

fl=  [0]

fr=  [0]

x=   [0]

abm= [0]

xl=  [0]

f_tilde = [0]

fex = [0]

i=1

for x in in_ff:

print i

a = x.split(" ")

for k in range(0,1):

print a[k]

xl=xl+[float(a[0])]

f=f+[float(a[1])]

fini=fini+[float(a[1])]

fn=fn+[float(a[1])]

fl=fl+[float(a[1])]

fr=fr+[float(a[1])]

f_tilde = f_tilde + [0]

fex = fex + [0]

df = df + [0]

dfn = dfn + [0]

i=i+1

imx=i

adm=[0.0]

for i in range(1,imx):

adm=adm+[0.0]

print "*"*40

print "xl,f"

for i in range(0,imx-1):

print float(xl[i]),"   ",float(f[i])

# set parameter ####################################

for n in range(0,1000):

###############################

for i in range(1,imx-1):

if u>=0.0:

ip = i - 1

dd = -dx

else:

ip = i + 1

dd = dx

xi = -u*dt

S  = (f[ip]-f[i])/dd

if S - f[ip] >0.0:

isgn = 1.0

else:

isgn = -1.0

TS = f[ip]-S

BB = ((abs(S-df[i])+eps)/(abs(S-df[ip])+eps)-1.0)/dd

a3 = (df[i]-S+(df[ip]-S)*(1.0+alpha*BB*dd))/(dd*dd)

a2 = S*alpha*BB+(S-df[i])/dd-a3*dd

a1 = df[i] +f[i]*alpha * BB

a0 = f[i]

fn[i]  = (a3*xi*xi*xi+a2*xi*xi+a1*xi+a0)/(1.0+alpha*BB*xi)

dfn[i] = (3.0*a3*xi*xi+2.0*a2*xi+a1)/(1.0+alpha*BB*xi)-(alpha*BB)*(fn[i])/(1.0+alpha*BB*xi)

#shift#########################################

for i in range(1,imx-1):

f[i]=fn[i]

df[i]=dfn[i]

print fn

###################################################

# Exact ###########################################

for i in range(0,imx-1):

if i>=50 and i<=60:

fex[i] = 1.0

else:

fex[i] = 0.0

#output############################################

ff=open("output.dat","w")

for i in range(0,imx-1):

ff.write(str(xl[i]))

ff.write("    ")

ff.write(str(f[i]))

ff.write("\n")

ff.close()

#output

#for i in range(0,imx-1):

# print f[i]

print "-"*100

print xl

print f

#####################graph#########################

y1=f

y2=fini

y3=fex

plot(xl,y1,label = 'RCIP')

plot(xl,y2,label = 'INITIAL')

plot(xl,y3,label = 'EXACT')

xlabel('x')

ylabel('f')

title('RCIP')

axis([-0,100.0,-0.2,1.2])

legend()

grid(True)

savefig('RCIP')

show()

###################################################