Skip to content

Commit 82b4401

Browse files
rbanerjee20rbanerjee20
committed
Added media folders for nonlin to fdm-jupyter-book.
1 parent da89af5 commit 82b4401

33 files changed

Lines changed: 3962 additions & 0 deletions

‎fdm-jupyter-book/notebooks/nonlin/exer-nonlin/fu_fem_int.py‎

Lines changed: 80 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,80 @@
1+
"""
2+
Explore algebraic forms arising from the integral f(u)*v in the finite
3+
element method.
4+
5+
| phi_im1 phi_i phi_ip1
6+
+1 /\ /\ /\
7+
| / \ / \ / \
8+
| / \ / \ / \
9+
| / \ / \ / \
10+
| / \ / \ / \
11+
| / \ / \ / \
12+
| / \ / \ / \
13+
| / \ / \ / \
14+
| / / \/ \
15+
| / / \ /\ \
16+
| / / \ / \ \
17+
| / / \ / \ \
18+
| / / \ / \ \
19+
| / / \ / \ \
20+
| / / \ / \ \
21+
| / / \ / \ \
22+
| / / \ / \ \
23+
--------------------------------------------------------------------------
24+
i-1 i i+1
25+
26+
cell L cell R
27+
"""
28+
from sympy import *
29+
import sys
30+
31+
x, u_im1, u_i, u_ip1, u, h, x_i = symbols('x u_im1 u_i u_ip1 u h x_i')
32+
33+
# Left cell: [x_im1, x_i]
34+
# Right cell: [x_i, x_ip1]
35+
x_im1 = x_i - h
36+
x_ip1 = x_i + h
37+
38+
phi = {'L': # Left cell
39+
{'im1': 1 - (x - x_im1)/h,
40+
'i': (x - x_im1)/h},
41+
'R': # Right cell
42+
{'i': 1 - (x-x_i)/h,
43+
'ip1': (x - x_i)/h}}
44+
45+
u = {'L': u_im1*phi['L']['im1'] + u_i *phi['L']['i'],
46+
'R': u_i *phi['R']['i'] + u_ip1*phi['R']['ip1']}
47+
48+
f = lambda u: eval(sys.argv[1])
49+
50+
integral_L = integrate(f(u['L'])*phi['L']['i'], (x, x_im1, x_i))
51+
integral_R = integrate(f(u['R'])*phi['R']['i'], (x, x_i, x_ip1))
52+
expr_i = simplify(expand(integral_L + integral_R))
53+
print expr_i
54+
latex_code = latex(expr_i, mode='plain')
55+
# Replace u_im1 sympy symbol name by latex symbol u_{i-1}
56+
latex_code = latex_code.replace('im1', '{i-1}')
57+
# Replace u_ip1 sympy symbol name by latex symbol u_{i+1}
58+
latex_code = latex_code.replace('ip1', '{i+1}')
59+
print latex_code
60+
# Escape (quote) latex_code so it can be sent as HTML text
61+
import cgi
62+
html_code = cgi.escape(latex_code)
63+
print html_code
64+
# Make a file with HTML code for displaying the LaTeX formula
65+
f = open('tmp.html', 'w')
66+
# Include an image that can be clicked on to yield a new
67+
# page with an interactive editor and display area where the
68+
# formula can be further edited
69+
text = """
70+
<a href="http://www.codecogs.com/eqnedit.php?latex=%(html_code)s"
71+
target="_blank">
72+
<img src="http://latex.codecogs.com/gif.latex?%(html_code)s"
73+
title="%(latex_code)s"/>
74+
</a>
75+
""" % vars()
76+
f.write(text)
77+
f.close()
78+
# load tmp.html into a browser
79+
80+

‎fdm-jupyter-book/notebooks/nonlin/exer-nonlin/logistic_p.py‎

Lines changed: 171 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,171 @@
1+
import numpy as np
2+
3+
def FE_logistic(p, u0, dt, Nt):
4+
u = np.zeros(Nt+1)
5+
u[0] = u0
6+
for n in range(Nt):
7+
u[n+1] = u[n] + dt*(1 - u[n])**p*u[n]
8+
return u
9+
10+
def BE_logistic(p, u0, dt, Nt, choice='Picard',
11+
eps_r=1E-3, omega=1, max_iter=1000):
12+
# u[n] = u[n-1] + dt*(1-u[n])**p*u[n]
13+
# -dt*(1-u[n])**p*u[n] + u[n] = u[n-1]
14+
if choice == 'Picard1':
15+
choice = 'Picard'
16+
max_iter = 1
17+
18+
u = np.zeros(Nt+1)
19+
iterations = []
20+
u[0] = u0
21+
for n in range(1, Nt+1):
22+
c = -u[n-1]
23+
if choice == 'Picard':
24+
def F(u):
25+
return -dt*(1-u)**p*u + u + c
26+
27+
u_ = u[n-1]
28+
k = 0
29+
while abs(F(u_)) > eps_r and k < max_iter:
30+
# u*(1-dt*(1-u_)**p) + c = 0
31+
u_ = omega*(-c/(1-dt*(1-u_)**p)) + (1-omega)*u_
32+
k += 1
33+
u[n] = u_
34+
iterations.append(k)
35+
36+
elif choice == 'Newton':
37+
def F(u):
38+
return -dt*(1-u)**p*u + u + c
39+
40+
def dF(u):
41+
return dt*p*(1-u)**(p-1)*u - dt*(1-u)**p + 1
42+
43+
u_ = u[n-1]
44+
k = 0
45+
while abs(F(u_)) > eps_r and k < max_iter:
46+
u_ = u_ - F(u_)/dF(u_)
47+
k += 1
48+
u[n] = u_
49+
iterations.append(k)
50+
return u, iterations
51+
52+
def CN_logistic(p, u0, dt, Nt):
53+
# u[n+1] = u[n] + dt*(1-u[n])**p*u[n+1]
54+
# (1 - dt*(1-u[n])**p)*u[n+1] = u[n]
55+
u = np.zeros(Nt+1)
56+
u[0] = u0
57+
for n in range(0, Nt):
58+
u[n+1] = u[n]/(1 - dt*(1 - u[n])**p)
59+
return u
60+
61+
def test_asymptotic_value():
62+
T = 100
63+
dt = 0.1
64+
Nt = int(round(T/float(dt)))
65+
u0 = 0.1
66+
p = 1.8
67+
68+
u_CN = CN_logistic(p, u0, dt, Nt)
69+
u_BE_Picard, iter_Picard = BE_logistic(
70+
p, u0, dt, Nt, choice='Picard',
71+
eps_r=1E-5, omega=1, max_iter=1000)
72+
u_BE_Newton, iter_Newton = BE_logistic(
73+
p, u0, dt, Nt, choice='Newton',
74+
eps_r=1E-5, omega=1, max_iter=1000)
75+
u_FE = FE_logistic(p, u0, dt, Nt)
76+
77+
for arr in u_CN, u_BE_Picard, u_BE_Newton, u_FE:
78+
expected = 1
79+
computed = arr[-1]
80+
tol = 0.01
81+
msg = 'expected=%s, computed=%s' % (expected, computed)
82+
print msg
83+
assert abs(expected - computed) < tol
84+
85+
from scitools.std import *
86+
87+
def demo():
88+
T = 12
89+
p = 1.2
90+
try:
91+
dt = float(sys.argv[1])
92+
eps_r = float(sys.argv[2])
93+
omega = float(sys.argv[3])
94+
except:
95+
dt = 0.8
96+
eps_r = 1E-3
97+
omega = 1
98+
N = int(round(T/float(dt)))
99+
100+
u_FE = FE_logistic(p, 0.1, dt, N)
101+
u_BE31, iter_BE31 = BE_logistic(p, 0.1, dt, N,
102+
'Picard1', eps_r, omega)
103+
u_BE3, iter_BE3 = BE_logistic(p, 0.1, dt, N,
104+
'Picard', eps_r, omega)
105+
u_BE4, iter_BE4 = BE_logistic(p, 0.1, dt, N,
106+
'Newton', eps_r, omega)
107+
u_CN = CN_logistic(p, 0.1, dt, N)
108+
109+
print 'Picard mean no of iterations (dt=%g):' % dt, \
110+
int(round(mean(iter_BE3)))
111+
print 'Newton mean no of iterations (dt=%g):' % dt, \
112+
int(round(mean(iter_BE4)))
113+
114+
t = np.linspace(0, dt*N, N+1)
115+
plot(t, u_FE, t, u_BE3, t, u_BE31, t, u_BE4, t, u_CN,
116+
legend=['FE', 'BE Picard', 'BE Picard1', 'BE Newton', 'CN gm'],
117+
title='dt=%g, eps=%.0E' % (dt, eps_r), xlabel='t', ylabel='u',
118+
legend_loc='lower right')
119+
filestem = 'logistic_N%d_eps%03d' % (N, log10(eps_r))
120+
savefig(filestem + '_u.png')
121+
savefig(filestem + '_u.pdf')
122+
figure()
123+
plot(range(1, len(iter_BE3)+1), iter_BE3, 'r-o',
124+
range(1, len(iter_BE4)+1), iter_BE4, 'b-o',
125+
legend=['Picard', 'Newton'],
126+
title='dt=%g, eps=%.0E' % (dt, eps_r),
127+
axis=[1, N+1, 0, max(iter_BE3 + iter_BE4)+1],
128+
xlabel='Time level', ylabel='No of iterations')
129+
savefig(filestem + '_iter.png')
130+
savefig(filestem + '_iter.pdf')
131+
raw_input()
132+
133+
def test_solvers():
134+
p = 2.5
135+
T = 5000
136+
dt = 0.5
137+
eps_r = 1E-6
138+
omega_values = [1]
139+
tol = 0.01
140+
N = int(round(T/float(dt)))
141+
142+
for omega in omega_values:
143+
u_FE = FE_logistic(p, 0.1, dt, N)
144+
u_BE31, iter_BE31 = BE_logistic(p, 0.1, dt, N,
145+
'Picard1', eps_r, omega)
146+
u_BE3, iter_BE3 = BE_logistic(p, 0.1, dt, N,
147+
'Picard', eps_r, omega)
148+
u_BE4, iter_BE4 = BE_logistic(p, 0.1, dt, N,
149+
'Newton', eps_r, omega)
150+
u_CN = CN_logistic(p, 0.1, dt, N)
151+
152+
print u_FE[-1], u_BE31[-1], u_BE3[-1], u_CN[-1]
153+
for u_x in u_FE, u_BE31, u_BE3, u_CN:
154+
print u_x[-1]
155+
assert abs(u_x[-1] - 1) < tol, 'u=%.16f' % u_x[-1]
156+
157+
"""
158+
t = np.linspace(0, dt*N, N+1)
159+
plot(t, u_FE, t, u_BE3, t, u_BE31, t, u_BE4, t, u_CN,
160+
legend=['FE', 'BE Picard', 'BE Picard1', 'BE Newton', 'CN gm'],
161+
title='dt=%g, eps=%.0E' % (dt, eps_r), xlabel='t', ylabel='u',
162+
legend_loc='lower right')
163+
filestem = 'tmp_N%d_eps%03d' % (N, log10(eps_r))
164+
savefig(filestem + '_u.png')
165+
savefig(filestem + '_u.pdf')
166+
"""
167+
168+
if __name__ == '__main__':
169+
#demo()
170+
#test_solvers()
171+
test_asymptotic_value()

‎fdm-jupyter-book/notebooks/nonlin/exer-nonlin/product_arith_mean_sympy.py‎

Lines changed: 33 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,33 @@
1+
from sympy import *
2+
t, dt = symbols('t dt')
3+
P, Q = symbols('P Q', cls=Function)
4+
5+
# Target expression P(t_{n+1/2})*Q(t_{n+1/2})
6+
# Simpler: P(0)*Q(0)
7+
# Arithmetic means of each factor:
8+
# 1/4*(P(-dt/2) + P(dt/2))*(Q(-dt/2) + Q(dt/2))
9+
# Arithmetic mean of the product:
10+
# 1/2*(P(-dt/2)*Q(-dt/2) + P(dt/2)*Q(dt/2))
11+
# Let's Taylor expand to compare
12+
13+
target = P(0)*Q(0)
14+
num_terms = 6
15+
P_p = P(t).series(t, 0, num_terms).subs(t, dt/2)
16+
print P_p
17+
P_m = P(t).series(t, 0, num_terms).subs(t, -dt/2)
18+
print P_m
19+
Q_p = Q(t).series(t, 0, num_terms).subs(t, dt/2)
20+
print Q_p
21+
Q_m = Q(t).series(t, 0, num_terms).subs(t, -dt/2)
22+
print Q_m
23+
24+
product_mean = Rational(1,2)*(P_m*Q_m + P_p*Q_p)
25+
product_mean = simplify(expand(product_mean))
26+
product_mean_error = product_mean - target
27+
28+
factor_mean = Rational(1,2)*(P_m + P_p)*Rational(1,2)*(Q_m + Q_p)
29+
factor_mean = simplify(expand(factor_mean))
30+
factor_mean_error = factor_mean - target
31+
32+
print 'product_mean_error:', product_mean_error
33+
print 'factor_mean_error:', factor_mean_error

‎fdm-jupyter-book/notebooks/nonlin/fig-nonlin/logistic_N100_eps-04_iter.pdf‎

14.1 KB
Binary file not shown.

‎fdm-jupyter-book/notebooks/nonlin/fig-nonlin/logistic_N100_eps-04_iter.png‎

24.9 KB
Loading

‎fdm-jupyter-book/notebooks/nonlin/fig-nonlin/logistic_N100_eps-04_u.pdf‎

22.8 KB
Binary file not shown.

‎fdm-jupyter-book/notebooks/nonlin/fig-nonlin/logistic_N100_eps-04_u.png‎

53.5 KB
Loading

‎fdm-jupyter-book/notebooks/nonlin/fig-nonlin/logistic_N10_eps-01_iter.pdf‎

13 KB
Binary file not shown.

‎fdm-jupyter-book/notebooks/nonlin/fig-nonlin/logistic_N10_eps-01_iter.png‎

32.4 KB
Loading

‎fdm-jupyter-book/notebooks/nonlin/fig-nonlin/logistic_N10_eps-01_u.pdf‎

16.9 KB
Binary file not shown.

0 commit comments

Comments
 (0)