print
print "RESULTS"
print
print
print "Figure 12"
print
var('y t')
mu = vector([2.4048255576957747,5.520078110286311,8.653727912911009,11.791534439014281,14.930917708487495,18.07106396791092,21.21163662987926,24.352471530748968, 27.493479132040253,30.634606468431976,33.77582021357357,36.91709835366399,40.05842576462824,43.19979171317677,46.34118837166139])
A = vector([0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0])
for i in range(0,15):
A[i]= 2*(mu[i] - N(pi)*struve_H(0,mu[i])/2.0)/mu[i]^2/bessel_J(1,mu[i])
print "KR = ", mu[0]^2, "(n = 1)"
print
def conc(x,t):
sum = 0.0
for j in range(0, 14):
sum = sum + A[j]*exp(-(mu[j]^2 - mu[0]^2)*t)*bessel_J(0,x*mu[j])
return sum
B = plot([],figsize=4,axes_labels=['Radial coordinate (r)','Concentration (c)'],gridlines=false, frame=True)
B += plot(x,(x,0,1),color='blue', legend_label='t = 0', ymin = 0, ymax = 1, linestyle = 'dashdot', thickness=1)
B += plot(conc(x,0.01),(x,0,1),color='green', legend_label='t = 0.01',ymin = 0, ymax = 1, linestyle = 'dotted', thickness=2)
B += plot(conc(x,0.05),(x,0,1),color='black', legend_label='t = 0.05', ymin = 0, ymax = 1, linestyle = 'dashed', thickness=1)
B += plot(conc(x,0.1),(x,0,1),color='red', legend_label='t = 0.1', ymin = 0, ymax = 1, linestyle = 'solid', thickness=1)
show(B)
print "KR = ", mu[1]^2, "(n = 2)"
print
def conc(x,t):
sum = 0.0
for j in range(1, 14):
sum = sum + A[j]*exp(-(mu[j]^2 - mu[1]^2)*t)*bessel_J(0,x*mu[j])
return sum
B = plot([],figsize=4,axes_labels=['Radial coordinate (r)','Concentration (c)'],gridlines=false, frame=True)
B += plot(x,(x,0,1),color='blue', legend_label='t = 0', ymin = 0, ymax = 1, linestyle = 'dashdot', thickness=1)
B += plot(conc(x,0.01),(x,0,1),color='green', legend_label='t = 0.01',ymin = 0, ymax = 1, linestyle = 'dotted', thickness=2)
B += plot(conc(x,0.05),(x,0,1),color='black', legend_label='t = 0.05', ymin = 0, ymax = 1, linestyle = 'dashed', thickness=1)
B += plot(conc(x,0.1),(x,0,1),color='red', legend_label='t = 0.1', ymin = 0, ymax = 1, linestyle = 'solid', thickness=1)
show(B)
print "KR = ", mu[2]^2, "(n = 3)"
print
def conc(x,t):
sum = 0.0
for j in range(2, 14):
sum = sum + A[j]*exp(-(mu[j]^2 - mu[2]^2)*t)*bessel_J(0,x*mu[j])
return sum
B = plot([],figsize=4,axes_labels=['Radial coordinate (r)','Concentration (c)'],gridlines=false, frame=True)
B += plot(x,(x,0,1),color='blue', legend_label='t = 0', ymin = 0, ymax = 1, linestyle = 'dashdot', thickness=1)
B += plot(conc(x,0.01),(x,0,1),color='green', legend_label='t = 0.01',ymin = 0, ymax = 1, linestyle = 'dotted', thickness=2)
B += plot(conc(x,0.05),(x,0,1),color='black', legend_label='t = 0.05', ymin = 0, ymax = 1, linestyle = 'dashed', thickness=1)
B += plot(conc(x,0.1),(x,0,1),color='red', legend_label='t = 0.1', ymin = 0, ymax = 1, linestyle = 'solid', thickness=1)
show(B)
print "KR = ", mu[3]^2, "(n = 4)"
print
def conc(x,t):
sum = 0.0
for j in range(3, 14):
sum = sum + A[j]*exp(-(mu[j]^2 - mu[3]^2)*t)*bessel_J(0,x*mu[j])
return sum
B = plot([],figsize=4,axes_labels=['Radial coordinate (r)','Concentration (c)'],gridlines=false, frame=True)
B += plot(x,(x,0,1),color='blue', legend_label='t = 0', ymin = 0, ymax = 1, linestyle = 'dashdot', thickness=1)
B += plot(conc(x,0.01),(x,0,1),color='green', legend_label='t = 0.01',ymin = 0, ymax = 1, linestyle = 'dotted', thickness=2)
B += plot(conc(x,0.05),(x,0,1),color='black', legend_label='t = 0.05', ymin = 0, ymax = 1, linestyle = 'dashed', thickness=1)
B += plot(conc(x,0.1),(x,0,1),color='red', legend_label='t = 0.1', ymin = 0, ymax = 1, linestyle = 'solid', thickness=1)
show(B)
print "KR = ", mu[4]^2, "(n = 5)"
print
def conc(x,t):
sum = 0.0
for j in range(4, 14):
sum = sum + A[j]*exp(-(mu[j]^2 - mu[4]^2)*t)*bessel_J(0,x*mu[j])
return sum
B = plot([],figsize=4,axes_labels=['Radial coordinate (r)','Concentration (c)'],gridlines=false, frame=True)
B += plot(x,(x,0,1),color='blue', legend_label='t = 0', ymin = 0, ymax = 1, linestyle = 'dashdot', thickness=1)
B += plot(conc(x,0.01),(x,0,1),color='green', legend_label='t = 0.01',ymin = 0, ymax = 1, linestyle = 'dotted', thickness=2)
B += plot(conc(x,0.05),(x,0,1),color='black', legend_label='t = 0.05', ymin = 0, ymax = 1, linestyle = 'dashed', thickness=1)
B += plot(conc(x,0.1),(x,0,1),color='red', legend_label='t = 0.1', ymin = 0, ymax = 1, linestyle = 'solid', thickness=1)
show(B)