Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

WIP: nonlin_ode #85

Open
wants to merge 4 commits into
base: master
Choose a base branch
from
Open

WIP: nonlin_ode #85

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
32a7df5
Change files to Python 3
jakubbober Sep 21, 2020
5ae98cb
General reformatting
jakubbober Sep 21, 2020
5dbdef8
Implement FE_logistic
jakubbober Sep 23, 2020
ad18d5f
Delete save parameter
jakubbober Sep 23, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
168 changes: 115 additions & 53 deletions fdm-devito-notebooks/05_nonlin/nonlin_ode.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -474,11 +474,23 @@
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [],
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\frac{dt - \\sqrt{dt^{2} + 4 dt u_{1} - 2 dt + 1} - 1}{2 dt}$"
],
"text/plain": [
"(dt - sqrt(dt**2 + 4*dt*u_1 - 2*dt + 1) - 1)/(2*dt)"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import sympy as sym\n",
"dt, u_1, u = sym.symbols('dt u_1 u')\n",
Expand All @@ -488,35 +500,59 @@
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [],
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\frac{dt + \\sqrt{dt^{2} + 4 dt u_{1} - 2 dt + 1} - 1}{2 dt}$"
],
"text/plain": [
"(dt + sqrt(dt**2 + 4*dt*u_1 - 2*dt + 1) - 1)/(2*dt)"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"r2"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [],
"execution_count": 18,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"-1/dt + 1 - u_1 + dt*(u_1**2 - u_1) + O(dt**2)\n"
]
}
],
"source": [
"print r1.series(dt, 0, 2) # 2 terms in dt, around dt=0"
"print(r1.series(dt, 0, 2)) # 2 terms in dt, around dt=0"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [],
"execution_count": 19,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"u_1 + dt*(-u_1**2 + u_1) + O(dt**2)\n"
]
}
],
"source": [
"print r2.series(dt, 0, 2)"
"print(r2.series(dt, 0, 2))"
]
},
{
Expand All @@ -540,24 +576,36 @@
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [],
"execution_count": 20,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"-oo\n"
]
}
],
"source": [
"print r1.limit(dt, 0)"
"print(r1.limit(dt, 0))"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false
},
"outputs": [],
"execution_count": 21,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"u_1\n"
]
}
],
"source": [
"print r2.limit(dt, 0)"
"print(r2.limit(dt, 0))"
]
},
{
Expand Down Expand Up @@ -624,11 +672,11 @@
"metadata": {},
"source": [
"Since the equation $\\hat F=0$ is only approximate, the solution $u$\n",
"does not equal the exact solution $\\uex$ of the exact\n",
"equation $F(\\uex)=0$, but we can hope that $u$ is closer to\n",
"$\\uex$ than $u^{-}$ is, and hence it makes sense to repeat the\n",
"does not equal the exact solution $u_\\text{e}$ of the exact\n",
"equation $F(u_\\text{e})=0$, but we can hope that $u$ is closer to\n",
"$u_\\text{e}$ than $u^{-}$ is, and hence it makes sense to repeat the\n",
"procedure, i.e., set $u^{-}=u$ and solve $\\hat F(u)=0$ again.\n",
"There is no guarantee that $u$ is closer to $\\uex$ than $u^{-}$,\n",
"There is no guarantee that $u$ is closer to $u_\\text{e}$ than $u^{-}$,\n",
"but this approach has proven to be effective in a wide range of\n",
"applications.\n",
"\n",
Expand Down Expand Up @@ -784,7 +832,7 @@
"\n",
"We shall later refer to the strategy of taking one Picard step, or\n",
"equivalently, linearizing terms with use of the solution at the\n",
"previous time step, as the *Picard1* method. mathcal{I}_t is a widely used\n",
"previous time step, as the *Picard1* method. It is a widely used\n",
"approach in science and technology, but with some limitations if\n",
"$\\Delta t$ is not sufficiently small (as will be illustrated later).\n",
"\n",
Expand Down Expand Up @@ -950,7 +998,7 @@
"The geometric mean approximation is often very effective for\n",
"linearizing quadratic nonlinearities. Both the arithmetic and geometric mean\n",
"approximations have truncation errors of order $\\Delta t^2$ and are\n",
"therefore compatible with the truncation error $\\Oof{\\Delta t^2}$\n",
"therefore compatible with the truncation error $\\mathcal{O}(\\Delta t^2)$\n",
"of the centered difference approximation for $u^\\prime$ in the Crank-Nicolson\n",
"method.\n",
"\n",
Expand Down Expand Up @@ -1066,7 +1114,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"mathcal{I}_t can be shown that the error in iteration $k+1$ of Newton's method is\n",
"It can be shown that the error in iteration $k+1$ of Newton's method is\n",
"proportional to\n",
"the square of the error in iteration $k$, a result referred to as\n",
"*quadratic convergence*. This means that for\n",
Expand All @@ -1076,7 +1124,7 @@
"However, the quadratic convergence appears only if $u^k$ is sufficiently\n",
"close to the solution. Further away from the solution the method can\n",
"easily converge very slowly or diverge. The reader is encouraged to do\n",
"[nonlin:exer:Newton:problems1](#nonlin:exer:Newton:problems1) to get a better understanding\n",
"[Problem 3](nonlin_exer.ipynb#nonlin:exer:Newton:problems1) to get a better understanding\n",
"for the behavior of the method.\n",
"\n",
"Application of Newton's method to the logistic equation discretized\n",
Expand Down Expand Up @@ -1229,7 +1277,7 @@
"<div id=\"nonlin:timediscrete:logistic:impl\"></div>\n",
"\n",
"\n",
"The program [`logistic.py`](${src_nonlin}/logistic.py) contains\n",
"The program [`logistic.py`](https://github.com/devitocodes/devito_book/blob/master/fdm-devito-notebooks/05_nonlin/src-nonlin/logistic.py) contains\n",
"implementations of all the methods described above.\n",
"Below is an extract of the file showing how the Picard and Newton\n",
"methods are implemented for a Backward Euler discretization of\n",
Expand All @@ -1241,9 +1289,7 @@
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false
},
"metadata": {},
"outputs": [],
"source": [
"def BE_logistic(u0, dt, Nt, choice='Picard',\n",
Expand Down Expand Up @@ -1302,9 +1348,7 @@
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": false
},
"metadata": {},
"outputs": [],
"source": [
"def CN_logistic(u0, dt, Nt):\n",
Expand Down Expand Up @@ -1336,7 +1380,7 @@
"Newton's method with one step comes far below the $\\epsilon_r$ tolerance,\n",
"while the Picard iteration needs on average 7 iterations to bring the\n",
"residual down to $\\epsilon_r=10^{-1}$, which gives insufficient\n",
"accuracy in the solution of the nonlinear equation. mathcal{I}_t is obvious\n",
"accuracy in the solution of the nonlinear equation. It is obvious\n",
"that the Picard1 method gives significant errors in addition to\n",
"the time discretization unless the time step is as small as in\n",
"the lower right plot.\n",
Expand Down Expand Up @@ -1549,7 +1593,7 @@
"but with $f(u,t)=\\sin(2(u+1))$, the $f(u^{-},t)u/u^{-}$ trick\n",
"leads to 7, 9, and 11 iterations during the first three steps, while\n",
"$f(u^{-},t)$ demands 17, 21, and 20 iterations.\n",
"(Experiments can be done with the code [`ODE_Picard_tricks.py`](${src_nonlin}/ODE_Picard_tricks.py).)\n",
"(Experiments can be done with the code [`ODE_Picard_tricks.py`](https://github.com/devitocodes/devito_book/blob/master/fdm-devito-notebooks/05_nonlin/src-nonlin/ODE_Picard_tricks.py)\n",
"\n",
"\n",
"\n",
Expand Down Expand Up @@ -2327,7 +2371,7 @@
"can be misleading when the initial start value for the iteration is\n",
"very close to the solution, since an unnecessary reduction in\n",
"the error measure is enforced. In such cases the absolute criteria\n",
"work better. mathcal{I}_t is common to combine the absolute and relative measures\n",
"work better. It is common to combine the absolute and relative measures\n",
"of the size of the residual,\n",
"as in"
]
Expand Down Expand Up @@ -2647,7 +2691,25 @@
]
}
],
"metadata": {},
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.2"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
6 changes: 3 additions & 3 deletions fdm-devito-notebooks/05_nonlin/src-nonlin/ODE_Picard_tricks.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,12 +24,12 @@ def f2(u, t):
tp = linspace(0, 4, 11)
u1, t = solver1.solve(tp)
u2, t = solver2.solve(tp)
print 'Picard it:', solver1.num_iterations_total
print 'Picard2 it:', solver2.num_iterations_total
print('Picard it:', solver1.num_iterations_total)
print('Picard2 it:', solver2.num_iterations_total)

from scitools.std import plot
plot(t, u1, t, u2, legend=('Picard', 'Picard2'))
raw_input()
input()

"""
f(u,t) = -u**3:
Expand Down
41 changes: 23 additions & 18 deletions fdm-devito-notebooks/05_nonlin/src-nonlin/logistic.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,11 +1,16 @@
import numpy as np

def FE_logistic(u0, dt, Nt):
u = np.zeros(N+1)
u[0] = u0
for n in range(Nt):
u[n+1] = u[n] + dt*(u[n] - u[n]**2)
return u
t = Dimension('t', spacing=Constant('h_t'))

u = TimeFunction(name='u', dimensions=(t,), shape=(Nt+1))
u.data[0] = u0

pde = u.dtr - u + u**2
stencil = Eq(u.forward, solve(pde, u.forward))
op = Operator(stencil)
op.apply(h_t=dt, t_M=Nt-1)
return u.data

def quadratic_roots(a, b, c):
delta = b**2 - 4*a*c
Expand Down Expand Up @@ -77,22 +82,22 @@ def quadratic_root_goes_to_infinity():
a = dt
b = 1 - dt
c = -0.1
print dt, quadratic_roots(a, b, c)
print(dt, quadratic_roots(a, b, c))

def sympy_analysis():
print 'sympy calculations'
print('sympy calculations')
import sympy as sym
dt, u_1, u = sym.symbols('dt u_1 u')
r1, r2 = sym.solve(dt*u**2 + (1-dt)*u - u_1, u)
print r1
print r2
print r1.series(dt, 0, 2)
print r2.series(dt, 0, 2)
print r1.limit(dt, 0)
print r2.limit(dt, 0)
print(r1)
print(r2)
print(r1.series(dt, 0, 2))
print(r2.series(dt, 0, 2))
print(r1.limit(dt, 0))
print(r2.limit(dt, 0))

sympy_analysis()
print '-----------------------------------------------------'
print('-----------------------------------------------------')
T = 9
try:
dt = float(sys.argv[1])
Expand All @@ -113,8 +118,8 @@ def sympy_analysis():
u_CN = CN_logistic(0.1, dt, N)

from numpy import mean
print 'Picard mean no of iterations (dt=%g):' % dt, int(round(mean(iter_BE3)))
print 'Newton mean no of iterations (dt=%g):' % dt, int(round(mean(iter_BE4)))
print('Picard mean no of iterations (dt=%g):' % dt, int(round(mean(iter_BE3))))
print('Newton mean no of iterations (dt=%g):' % dt, int(round(mean(iter_BE4))))

t = np.linspace(0, dt*N, N+1)
plot(t, u_FE, t, u_BE2, t, u_BE3, t, u_BE31, t, u_BE4, t, u_CN,
Expand All @@ -125,8 +130,8 @@ def sympy_analysis():
savefig(filestem + '_u.png')
savefig(filestem + '_u.pdf')
figure()
plot(range(1, len(iter_BE3)+1), iter_BE3, 'r-o',
range(1, len(iter_BE4)+1), iter_BE4, 'b-o',
plot(list(range(1, len(iter_BE3)+1)), iter_BE3, 'r-o',
list(range(1, len(iter_BE4)+1)), iter_BE4, 'b-o',
legend=['Picard', 'Newton'], title='dt=%g, eps=%.0E' % (dt, eps_r),
axis=[1, N+1, 0, max(iter_BE3 + iter_BE4)+1],
xlabel='Time level', ylabel='No of iterations')
Expand Down
Loading