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()
###################################################
