{
"cells": [
{
"cell_type": "markdown",
"id": "407bf032",
"metadata": {},
"source": [
"# Error Estimates and Stability"
]
},
{
"cell_type": "markdown",
"id": "29147f32",
"metadata": {},
"source": [
"This notebook introduces the estimation of errors via the combination of quadrature rules. The second section illustrates the stability of Euler's method. In the third section, the backward Euler method is illustrated on an example of a stiff equation."
]
},
{
"cell_type": "markdown",
"id": "8a0332c8",
"metadata": {},
"source": [
"The last section illustrates the application of the ode solver with the DifferentialEquations package."
]
},
{
"cell_type": "markdown",
"id": "8832582c",
"metadata": {},
"source": [
"## 1. Estimating the Error of the Trapezoidal Rule"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "ade3f6c0",
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"\n",
"(process:61192): GLib-GIO-WARNING **: 19:04:39.189: Unexpectedly, UWP app `Clipchamp.Clipchamp_2.5.1.0_neutral__yxz26nhyzhsrt' (AUMId `Clipchamp.Clipchamp_yxz26nhyzhsrt!App') supports 41 extensions but has no verbs\n",
"\n",
"(process:61192): GLib-GIO-WARNING **: 19:04:39.220: Unexpectedly, UWP app `Evernote.Evernote_10.47.7.0_x64__q4d96b2w5wcc2' (AUMId `Evernote.Evernote_q4d96b2w5wcc2!Evernote') supports 1 extensions but has no verbs\n"
]
}
],
"source": [
"using Plots"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "194e197b",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"f (generic function with 1 method)"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f(x) = (x-1)^3 + 0.1*x^2 - x + 2"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "63c17404",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"r = -0.35:0.01:2.2\n",
"y = [f(t) for t in r]\n",
"plot(r,y,label=\"y=f(x)\",legend=:topleft)"
]
},
{
"cell_type": "markdown",
"id": "74e2c66f",
"metadata": {},
"source": [
"## The Trapezoidal Rule"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "d2ed66e8",
"metadata": {},
"outputs": [],
"source": [
"plot(r,y,label=\"y=f(x)\",legend=:topleft)\n",
"intrange = 0:0.01:2\n",
"mf0 = [0 for t in intrange]\n",
"plot!(intrange, mf0, line=(:red), label=\"\")\n",
"yT(x) = f(0) + (f(2)-f(0))*x/2\n",
"my = [yT(t) for t in intrange]\n",
"plot!(intrange, my, line=(:red), label=\"\")\n",
"plot!([0.0; 0.0], [0.0; f(0)], line=(:red), label=\"\")\n",
"plot!([2.0; 2.0], [0.0; f(2)], line=(:red), label=\"\")\n",
"plot!([0.0; 0.0], [0.0; f(0)], line=(:black, :dot), label=\"\")\n",
"plot!([2.0; 2.0], [0.0; f(2)], line=(:black, :dot), label=\"\")\n",
"scatter!([0.0; 0.0],[f(0); f(0)],marker=(:blue),label=\"f(0)\")\n",
"scatter!([2.0; 2.0],[f(2); f(2)],marker=(:green),label=\"f(2)\")\n",
"savefig(\"figtrapezoidalrule.png\")"
]
},
{
"cell_type": "markdown",
"id": "60652232",
"metadata": {},
"source": [
"## Simpson's Rule"
]
},
{
"cell_type": "markdown",
"id": "6cd1844c",
"metadata": {},
"source": [
"We construct Simpson's rule on the interval $[0,2]$, evaluating ``f`` at 0, 1 and 2. The weights are the integrals of the lagrange polynomials of the quadratic polynomial interpolating at 0, 1, and 2."
]
},
{
"cell_type": "markdown",
"id": "119edb5b",
"metadata": {},
"source": [
"To compute the weights, we first construct the Lagrange polynomials, using ``SymPy``."
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "180da52f",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$x$"
],
"text/plain": [
"x"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using SymPy\n",
"x = Sym(\"x\")"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "f2508f73",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The first Lagrange polynomial : (x - 2)*(x - 1)/2\n"
]
},
{
"data": {
"text/latex": [
"$\\frac{1}{3}$"
],
"text/plain": [
"1/3"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"L0 = (x-1)*(x-2)/((0-1)*(0-2))\n",
"println(\"The first Lagrange polynomial : \", L0)\n",
"w0 = integrate(L0, (x, 0, 2))"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "4976ed40",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The second Lagrange polynomial : -x*(x - 2)\n"
]
},
{
"data": {
"text/latex": [
"$\\frac{4}{3}$"
],
"text/plain": [
"4/3"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"L1 = (x-0)*(x-2)/((1-0)*(1-2))\n",
"println(\"The second Lagrange polynomial : \", L1)\n",
"w1 = integrate(L1, (x, 0, 2))"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "71b367f1",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The third Lagrange polynomial : x*(x - 1)/2\n"
]
},
{
"data": {
"text/latex": [
"$\\frac{1}{3}$"
],
"text/plain": [
"1/3"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"L2 = (x-0)*(x-1)/((2-0)*(2-1))\n",
"println(\"The third Lagrange polynomial : \", L2)\n",
"w2 = integrate(L2, (x, 0, 2))"
]
},
{
"cell_type": "markdown",
"id": "7e1863c1",
"metadata": {},
"source": [
"The output of the three cells above are the three weights of Simpson's rule. With the Lagrange polynomials we also have an expression for the interpolating polynomial ``p``, interpolating at ``(0, f(0))``, ``(1, f(1))``, and ``(2, f(2))``."
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "6143ee70",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$- 1.1 x \\left(x - 2\\right) + 0.7 x \\left(x - 1\\right) + 0.5 \\left(x - 2\\right) \\left(x - 1\\right)$"
],
"text/plain": [
"-1.1*x*(x - 2) + 0.7*x*(x - 1) + 0.5*(x - 2)*(x - 1)"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"p = f(0)*L0 + f(1)*L1 + f(2)*L2"
]
},
{
"cell_type": "markdown",
"id": "c3de09c6",
"metadata": {},
"source": [
"To verify the interpolation conditions and to make the plot, we turn the expression for ``p`` into a Julia function."
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "dad7d2e9",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"#118 (generic function with 1 method)"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"fp = lambdify(p)"
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "baf7d8a5",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"fp(0) : 1.0 f(0) : 1.0\n",
"fp(1) : 1.1 f(1) : 1.1\n",
"fp(2) : 1.4 f(2) : 1.4\n"
]
}
],
"source": [
"println(\"fp(0) : \", fp(0), \" f(0) : \", f(0))\n",
"println(\"fp(1) : \", fp(1), \" f(1) : \", f(1))\n",
"println(\"fp(2) : \", fp(2), \" f(2) : \", f(2))"
]
},
{
"cell_type": "markdown",
"id": "be4b493c",
"metadata": {},
"source": [
"Now we can make the plot of Simpson's rule."
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "edd12073",
"metadata": {},
"outputs": [],
"source": [
"plot(r,y,label=\"y=f(x)\",legend=:topleft)\n",
"intrange = 0:0.01:2\n",
"mf0 = [0 for t in intrange]\n",
"plot!(intrange, mf0, line=(:green), label=\"\")\n",
"my = [fp(t) for t in intrange]\n",
"plot!(intrange, my, line=(:green), label=\"\")\n",
"plot!([0.0; 0.0], [0.0; f(0)], line=(:green), label=\"\")\n",
"plot!([2.0; 2.0], [0.0; f(2)], line=(:green), label=\"\")\n",
"plot!([0.0; 0.0], [0.0; f(0)], line=(:black, :dot), label=\"\")\n",
"plot!([1.0; 1.0], [0.0; f(1)], line=(:black, :dot), label=\"\")\n",
"plot!([2.0; 2.0], [0.0; f(2)], line=(:black, :dot), label=\"\")\n",
"scatter!([0.0; 0.0],[f(0); f(0)],marker=(:blue),label=\"f(0)\")\n",
"scatter!([1.0; 1.0],[f(1); f(1)],marker=(:purple),label=\"f(1)\")\n",
"scatter!([2.0; 2.0],[f(2); f(2)],marker=(:green),label=\"f(2)\")\n",
"savefig(\"figsimpsonsrule.png\")"
]
},
{
"cell_type": "markdown",
"id": "ea5133c4",
"metadata": {},
"source": [
"Now we can visualize the idea of the estimation of the error of the Trapezoidal rule. Simpson's rule uses one extra function evaluation in the middle of the interval. Comparing the value of Simpson's rule with the value of the Trapezoidal rule gives an estimate for the error of the Trapezoidal rule, at the expense of one extra function evaluation."
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "b15b7a4b",
"metadata": {},
"outputs": [],
"source": [
"plot(r,y,label=\"y=f(x)\",legend=:topleft)\n",
"intrange = 0:0.01:2\n",
"mf0 = [0 for t in intrange]\n",
"plot!(intrange, mf0, line=(:red), label=\"\")\n",
"yT(x) = f(0) + (f(2)-f(0))*x/2\n",
"my = [yT(t) for t in intrange]\n",
"plot!(intrange, my, line=(:red), label=\"\")\n",
"plot!([0.0; 0.0], [0.0; f(0)], line=(:red), label=\"\")\n",
"plot!([2.0; 2.0], [0.0; f(2)], line=(:red), label=\"\")\n",
"plot!([0.0; 0.0], [0.0; f(0)], line=(:black, :dot), label=\"\")\n",
"# augment with the plot for Simpson's rule\n",
"mf0 = [0 for t in intrange]\n",
"plot!(intrange, mf0, line=(:green), label=\"\")\n",
"my = [fp(t) for t in intrange]\n",
"plot!(intrange, my, line=(:green), label=\"\")\n",
"plot!([0.0; 0.0], [0.0; f(0)], line=(:green), label=\"\")\n",
"plot!([2.0; 2.0], [0.0; f(2)], line=(:green), label=\"\")\n",
"plot!([0.0; 0.0], [0.0; f(0)], line=(:black, :dot), label=\"\")\n",
"plot!([1.0; 1.0], [0.0; f(1)], line=(:black, :dot), label=\"\")\n",
"plot!([2.0; 2.0], [0.0; f(2)], line=(:black, :dot), label=\"\")\n",
"plot!([2.0; 2.0], [0.0; f(2)], line=(:black, :dot), label=\"\")\n",
"scatter!([0.0; 0.0],[f(0); f(0)],marker=(:blue),label=\"f(0)\")\n",
"scatter!([1.0; 1.0],[f(1); f(1)],marker=(:purple),label=\"f(1)\")\n",
"scatter!([2.0; 2.0],[f(2); f(2)],marker=(:green),label=\"f(2)\")\n",
"savefig(\"figtrapeziumsimpson.png\")"
]
},
{
"cell_type": "markdown",
"id": "5b472960",
"metadata": {},
"source": [
"## 2. The Stability of Euler's method"
]
},
{
"cell_type": "markdown",
"id": "29f9d355",
"metadata": {},
"source": [
"Let us examine the behavior of Euler's method on problems where the solution is exponentially decaying."
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "4e127e20",
"metadata": {},
"outputs": [],
"source": [
"using Printf"
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "8bde9f8e",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"eulerexp"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"\"\"\"\n",
" eulerexp(n::Int64,lambda::Float64,verbose::Bool=true)\n",
"\n",
"Euler's method with n evaluations\n",
"on the interval [0,1] on the test\n",
"equation y' = lambda*y, y(0) = 1.\n",
"\n",
"On return is a tuple with two arrays, each with n+1 values,\n",
"for the approximations y at the corresponding x values.\n",
"\"\"\"\n",
"function eulerexp(n::Int64,lambda::Float64,verbose::Bool=true)\n",
" h = 1.0/n\n",
" if verbose\n",
" println(\" i x approx exact error\")\n",
" end\n",
" xv = zeros(n+1)\n",
" yv = zeros(n+1)\n",
" y = 0.0\n",
" for i=0:n\n",
" x = i*h\n",
" y = (1.0+lambda*h)^i\n",
" (xv[i+1], yv[i+1]) = (x, y)\n",
" exact = exp(lambda*x)\n",
" if verbose\n",
" stri = @sprintf(\"%3d\", i)\n",
" strx = @sprintf(\"%.2f\", x)\n",
" stry = @sprintf(\"%.6e\", y)\n",
" strexp = @sprintf(\"%.6e\", exact)\n",
" strerr = @sprintf(\"%.2e\", abs(y-exact))\n",
" println(\"$stri $strx $stry $strexp $strerr\")\n",
" end\n",
" end\n",
" return (xv, yv)\n",
"end"
]
},
{
"cell_type": "markdown",
"id": "a7444381",
"metadata": {},
"source": [
"For testing purposes, we do a run with ``lambda = 1.0``."
]
},
{
"cell_type": "code",
"execution_count": 16,
"id": "775b713f",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" i x approx exact error\n",
" 0 0.00 1.000000e+00 1.000000e+00 0.00e+00\n",
" 1 0.10 1.100000e+00 1.105171e+00 5.17e-03\n",
" 2 0.20 1.210000e+00 1.221403e+00 1.14e-02\n",
" 3 0.30 1.331000e+00 1.349859e+00 1.89e-02\n",
" 4 0.40 1.464100e+00 1.491825e+00 2.77e-02\n",
" 5 0.50 1.610510e+00 1.648721e+00 3.82e-02\n",
" 6 0.60 1.771561e+00 1.822119e+00 5.06e-02\n",
" 7 0.70 1.948717e+00 2.013753e+00 6.50e-02\n",
" 8 0.80 2.143589e+00 2.225541e+00 8.20e-02\n",
" 9 0.90 2.357948e+00 2.459603e+00 1.02e-01\n",
" 10 1.00 2.593742e+00 2.718282e+00 1.25e-01\n"
]
},
{
"data": {
"text/plain": [
"([0.0, 0.1, 0.2, 0.30000000000000004, 0.4, 0.5, 0.6000000000000001, 0.7000000000000001, 0.8, 0.9, 1.0], [1.0, 1.1, 1.2100000000000002, 1.3310000000000004, 1.4641000000000004, 1.6105100000000006, 1.7715610000000008, 1.9487171000000012, 2.1435888100000016, 2.357947691000002, 2.5937424601000023])"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"(x, y) = eulerexp(10, 1.0)"
]
},
{
"cell_type": "code",
"execution_count": 17,
"id": "02a500af",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"xr = 0.0:0.01:1.0\n",
"yr = [exp(t) for t in xr]\n",
"plot(xr,yr, label=\"y = exp(x)\", legend=:topleft)\n",
"scatter!(x,y, label=\"euler\")"
]
},
{
"cell_type": "markdown",
"id": "ab1ab4c0",
"metadata": {},
"source": [
"The plot is okay. Now let us consider an exponentially decaying function, for ``lambda = -40``."
]
},
{
"cell_type": "code",
"execution_count": 18,
"id": "5733f8a0",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" i x approx exact error\n",
" 0 0.00 1.000000e+00 1.000000e+00 0.00e+00\n",
" 1 0.10 -3.000000e+00 1.831564e-02 3.02e+00\n",
" 2 0.20 9.000000e+00 3.354626e-04 9.00e+00\n",
" 3 0.30 -2.700000e+01 6.144212e-06 2.70e+01\n",
" 4 0.40 8.100000e+01 1.125352e-07 8.10e+01\n",
" 5 0.50 -2.430000e+02 2.061154e-09 2.43e+02\n",
" 6 0.60 7.290000e+02 3.775135e-11 7.29e+02\n",
" 7 0.70 -2.187000e+03 6.914400e-13 2.19e+03\n",
" 8 0.80 6.561000e+03 1.266417e-14 6.56e+03\n",
" 9 0.90 -1.968300e+04 2.319523e-16 1.97e+04\n",
" 10 1.00 5.904900e+04 4.248354e-18 5.90e+04\n"
]
}
],
"source": [
"(x,y) = eulerexp(10,-40.0);"
]
},
{
"cell_type": "markdown",
"id": "78f9aed2",
"metadata": {},
"source": [
"We see that the solution diverges!"
]
},
{
"cell_type": "code",
"execution_count": 19,
"id": "5060ab81",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" i x approx exact error\n",
" 0 0.00 1.000000e+00 1.000000e+00 0.00e+00\n",
" 1 0.05 -1.000000e+00 1.353353e-01 1.14e+00\n",
" 2 0.10 1.000000e+00 1.831564e-02 9.82e-01\n",
" 3 0.15 -1.000000e+00 2.478752e-03 1.00e+00\n",
" 4 0.20 1.000000e+00 3.354626e-04 1.00e+00\n",
" 5 0.25 -1.000000e+00 4.539993e-05 1.00e+00\n",
" 6 0.30 1.000000e+00 6.144212e-06 1.00e+00\n",
" 7 0.35 -1.000000e+00 8.315287e-07 1.00e+00\n",
" 8 0.40 1.000000e+00 1.125352e-07 1.00e+00\n",
" 9 0.45 -1.000000e+00 1.522998e-08 1.00e+00\n",
" 10 0.50 1.000000e+00 2.061154e-09 1.00e+00\n",
" 11 0.55 -1.000000e+00 2.789468e-10 1.00e+00\n",
" 12 0.60 1.000000e+00 3.775135e-11 1.00e+00\n",
" 13 0.65 -1.000000e+00 5.109089e-12 1.00e+00\n",
" 14 0.70 1.000000e+00 6.914400e-13 1.00e+00\n",
" 15 0.75 -1.000000e+00 9.357623e-14 1.00e+00\n",
" 16 0.80 1.000000e+00 1.266417e-14 1.00e+00\n",
" 17 0.85 -1.000000e+00 1.713908e-15 1.00e+00\n",
" 18 0.90 1.000000e+00 2.319523e-16 1.00e+00\n",
" 19 0.95 -1.000000e+00 3.139133e-17 1.00e+00\n",
" 20 1.00 1.000000e+00 4.248354e-18 1.00e+00\n"
]
}
],
"source": [
"(x,y) = eulerexp(20,-40.0);\n",
"xr = 0.0:0.01:1.0\n",
"yr = [exp(-40*t) for t in xr]\n",
"plot(xr,yr, label=\"y = exp(-40x)\", legend=3)\n",
"scatter!(x,y, label=\"euler\")\n",
"savefig(\"figeuleroscillates.png\")"
]
},
{
"cell_type": "markdown",
"id": "fe279a4c",
"metadata": {},
"source": [
"## 3. An Example of a Stiff Equation"
]
},
{
"cell_type": "markdown",
"id": "ea231f35",
"metadata": {},
"source": [
"The equation $y' = 10(1-y)$, $y(0) = 0.5$ is an example of a stiff equation, with exact solution $y(x) = 1 - \\exp(-10x)/2$."
]
},
{
"cell_type": "code",
"execution_count": 20,
"id": "bcd0728a",
"metadata": {},
"outputs": [],
"source": [
"r = 0.0:0.01:1.0\n",
"yf(x) = 1.0 - exp(-10*x)/2\n",
"yr = [yf(t) for t in r]\n",
"plot(r,yr,label=\"1-exp(10x)/2\",legend=3)\n",
"savefig(\"figstiffequation.png\")"
]
},
{
"cell_type": "code",
"execution_count": 21,
"id": "c0cf7891",
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"WARNING: using SciPy.integrate in module Main conflicts with an existing identifier.\n"
]
}
],
"source": [
"using SciPy"
]
},
{
"cell_type": "markdown",
"id": "6b93cdf0",
"metadata": {},
"source": [
"The ``integrate`` of SciPy conflicts with the ``integrate`` of SymPy."
]
},
{
"cell_type": "code",
"execution_count": 22,
"id": "75d07639",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"PyObject "
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"\"\"\"\n",
" rhs(t, y)\n",
"\n",
"Defines the right hand side of the ODE.\n",
"\"\"\"\n",
"function rhs(t, y)\n",
" result = zeros(length(y))\n",
" result[1] = 10*(1.0 - y[1])\n",
" return result\n",
"end\n",
"\n",
"sol = SciPy.integrate.RK45(rhs, 0.0, (0.5, ), 100.0, rtol=1.0e-8)"
]
},
{
"cell_type": "code",
"execution_count": 23,
"id": "7e120f64",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"step : h : t : y : error\n",
" 1 : 1.15e-02 : 1.15e-02 : 5.5431000898772287e-01 : 3.47e-10\n",
" 2 : 3.84e-02 : 2.69e-02 : 6.5936065774270014e-01 : 5.56e-08\n",
" 3 : 6.56e-02 : 2.72e-02 : 7.4043603599393226e-01 : 8.77e-08\n",
" 4 : 9.42e-02 : 2.87e-02 : 8.0515781530443864e-01 : 1.14e-07\n",
" 5 : 1.24e-01 : 3.03e-02 : 8.5602147176151033e-01 : 1.34e-07\n",
" 6 : 1.56e-01 : 3.20e-02 : 8.9545342045424714e-01 : 1.50e-07\n",
" 7 : 1.90e-01 : 3.40e-02 : 9.2555426010467901e-01 : 1.62e-07\n",
" 8 : 2.27e-01 : 3.61e-02 : 9.4813813668934199e-01 : 1.70e-07\n",
" 9 : 2.65e-01 : 3.86e-02 : 9.6475494222001734e-01 : 1.76e-07\n",
" 10 : 3.07e-01 : 4.14e-02 : 9.7671343059300542e-01 : 1.80e-07\n",
" 11 : 3.51e-01 : 4.47e-02 : 9.8510404980970256e-01 : 1.82e-07\n",
" 12 : 4.00e-01 : 4.84e-02 : 9.9082145912026043e-01 : 1.83e-07\n",
" 13 : 4.53e-01 : 5.28e-02 : 9.9458662631764372e-01 : 1.83e-07\n",
" 14 : 5.11e-01 : 5.80e-02 : 9.9696839982757801e-01 : 1.81e-07\n",
" 15 : 5.75e-01 : 6.42e-02 : 9.9840443896301689e-01 : 1.79e-07\n",
" 16 : 6.47e-01 : 7.18e-02 : 9.9922136841700548e-01 : 1.76e-07\n",
" 17 : 7.28e-01 : 8.11e-02 : 9.9965399850011893e-01 : 1.73e-07\n",
" 18 : 8.21e-01 : 9.30e-02 : 9.9986342001049755e-01 : 1.67e-07\n",
" 19 : 9.29e-01 : 1.08e-01 : 9.9995374183362096e-01 : 1.61e-07\n",
" 20 : 1.06e+00 : 1.29e-01 : 9.9998719242972067e-01 : 1.52e-07\n",
" 21 : 1.22e+00 : 1.58e-01 : 9.9999726064227368e-01 : 1.40e-07\n",
" 22 : 1.42e+00 : 2.01e-01 : 9.9999952542410586e-01 : 1.25e-07\n",
" 23 : 1.68e+00 : 2.66e-01 : 9.9999985355481424e-01 : 1.22e-07\n",
" 24 : 2.05e+00 : 3.64e-01 : 9.9999973485715155e-01 : 2.65e-07\n",
" 25 : 2.42e+00 : 3.75e-01 : 9.9999941869825737e-01 : 5.81e-07\n",
" 26 : 2.75e+00 : 3.30e-01 : 9.9999942364097416e-01 : 5.76e-07\n",
" 27 : 3.08e+00 : 3.30e-01 : 9.9999942854166390e-01 : 5.71e-07\n",
" 28 : 3.42e+00 : 3.38e-01 : 9.9999934310534966e-01 : 6.57e-07\n",
" 29 : 3.76e+00 : 3.37e-01 : 9.9999925608084816e-01 : 7.44e-07\n",
" 30 : 4.09e+00 : 3.28e-01 : 9.9999928729551646e-01 : 7.13e-07\n",
" 31 : 4.41e+00 : 3.22e-01 : 9.9999939532208404e-01 : 6.05e-07\n",
" 32 : 4.73e+00 : 3.26e-01 : 9.9999944878319336e-01 : 5.51e-07\n",
" 33 : 5.07e+00 : 3.36e-01 : 9.9999939390138559e-01 : 6.06e-07\n",
" 34 : 5.41e+00 : 3.40e-01 : 9.9999927698887037e-01 : 7.23e-07\n",
" 35 : 5.74e+00 : 3.33e-01 : 9.9999924290782694e-01 : 7.57e-07\n",
" 36 : 6.07e+00 : 3.23e-01 : 9.9999934518710076e-01 : 6.55e-07\n",
" 37 : 6.39e+00 : 3.22e-01 : 9.9999944660266149e-01 : 5.53e-07\n",
" 38 : 6.72e+00 : 3.31e-01 : 9.9999944068942237e-01 : 5.59e-07\n",
" 39 : 7.06e+00 : 3.41e-01 : 9.9999932632415245e-01 : 6.74e-07\n",
" 40 : 7.40e+00 : 3.38e-01 : 9.9999922296192101e-01 : 7.77e-07\n",
" 41 : 7.73e+00 : 3.27e-01 : 9.9999928015266826e-01 : 7.20e-07\n",
" 42 : 8.04e+00 : 3.19e-01 : 9.9999941691195637e-01 : 5.83e-07\n",
" 43 : 8.37e+00 : 3.25e-01 : 9.9999947069777839e-01 : 5.29e-07\n",
" 44 : 8.71e+00 : 3.38e-01 : 9.9999939102462720e-01 : 6.09e-07\n",
" 45 : 9.05e+00 : 3.43e-01 : 9.9999924225315584e-01 : 7.58e-07\n",
" 46 : 9.38e+00 : 3.32e-01 : 9.9999921677931169e-01 : 7.83e-07\n",
" 47 : 9.70e+00 : 3.20e-01 : 9.9999935811293950e-01 : 6.42e-07\n",
" 48 : 1.00e+01 : 3.20e-01 : 9.9999947434112868e-01 : 5.26e-07\n",
" 49 : 1.04e+01 : 3.33e-01 : 9.9999945274654811e-01 : 5.47e-07\n",
" 50 : 1.07e+01 : 3.44e-01 : 9.9999930123715197e-01 : 6.99e-07\n",
" 51 : 1.10e+01 : 3.39e-01 : 9.9999918044910197e-01 : 8.20e-07\n",
" 52 : 1.14e+01 : 3.24e-01 : 9.9999927549895740e-01 : 7.25e-07\n",
" 53 : 1.17e+01 : 3.16e-01 : 9.9999944474678792e-01 : 5.55e-07\n",
" 54 : 1.20e+01 : 3.25e-01 : 9.9999949517552722e-01 : 5.05e-07\n",
" 55 : 1.23e+01 : 3.42e-01 : 9.9999938331547633e-01 : 6.17e-07\n",
" 56 : 1.27e+01 : 3.45e-01 : 9.9999919481762745e-01 : 8.05e-07\n",
" 57 : 1.30e+01 : 3.31e-01 : 9.9999918712252944e-01 : 8.13e-07\n",
" 58 : 1.33e+01 : 3.16e-01 : 9.9999937765194880e-01 : 6.22e-07\n",
" 59 : 1.37e+01 : 3.18e-01 : 9.9999950680793137e-01 : 4.93e-07\n",
" 60 : 1.40e+01 : 3.35e-01 : 9.9999946379699600e-01 : 5.36e-07\n",
" 61 : 1.43e+01 : 3.48e-01 : 9.9999926462238686e-01 : 7.35e-07\n",
" 62 : 1.47e+01 : 3.40e-01 : 9.9999912647460543e-01 : 8.74e-07\n",
" 63 : 1.50e+01 : 3.21e-01 : 9.9999927510839381e-01 : 7.25e-07\n",
" 64 : 1.53e+01 : 3.13e-01 : 9.9999947918956966e-01 : 5.21e-07\n",
" 65 : 1.56e+01 : 3.26e-01 : 9.9999952161088235e-01 : 4.78e-07\n",
" 66 : 1.60e+01 : 3.46e-01 : 9.9999936859160299e-01 : 6.31e-07\n",
" 67 : 1.63e+01 : 3.48e-01 : 9.9999913055419620e-01 : 8.69e-07\n",
" 68 : 1.67e+01 : 3.29e-01 : 9.9999915475742018e-01 : 8.45e-07\n",
" 69 : 1.70e+01 : 3.12e-01 : 9.9999940496327122e-01 : 5.95e-07\n",
" 70 : 1.73e+01 : 3.16e-01 : 9.9999954342845954e-01 : 4.57e-07\n",
" 71 : 1.76e+01 : 3.39e-01 : 9.9999947262655153e-01 : 5.27e-07\n",
" 72 : 1.80e+01 : 3.53e-01 : 9.9999921222585875e-01 : 7.88e-07\n",
" 73 : 1.83e+01 : 3.40e-01 : 9.9999905900672814e-01 : 9.41e-07\n",
" 74 : 1.86e+01 : 3.16e-01 : 9.9999928101099600e-01 : 7.19e-07\n",
" 75 : 1.89e+01 : 3.09e-01 : 9.9999951985441415e-01 : 4.80e-07\n",
" 76 : 1.93e+01 : 3.27e-01 : 9.9999954903434196e-01 : 4.51e-07\n",
" 77 : 1.96e+01 : 3.51e-01 : 9.9999934408598190e-01 : 6.56e-07\n",
" 78 : 2.00e+01 : 3.51e-01 : 9.9999904444622090e-01 : 9.56e-07\n",
" 79 : 2.03e+01 : 3.26e-01 : 9.9999912132968738e-01 : 8.79e-07\n",
" 80 : 2.06e+01 : 3.06e-01 : 9.9999944044591815e-01 : 5.60e-07\n",
" 81 : 2.09e+01 : 3.15e-01 : 9.9999958295674396e-01 : 4.17e-07\n",
" 82 : 2.13e+01 : 3.43e-01 : 9.9999947777162757e-01 : 5.22e-07\n",
" 83 : 2.16e+01 : 3.59e-01 : 9.9999913875351398e-01 : 8.61e-07\n",
" 84 : 2.20e+01 : 3.40e-01 : 9.9999897661548121e-01 : 1.02e-06\n",
" 85 : 2.23e+01 : 3.11e-01 : 9.9999929501137819e-01 : 7.05e-07\n",
" 86 : 2.26e+01 : 3.05e-01 : 9.9999956528214817e-01 : 4.35e-07\n",
" 87 : 2.29e+01 : 3.29e-01 : 9.9999957613682422e-01 : 4.24e-07\n",
" 88 : 2.33e+01 : 3.58e-01 : 9.9999930659970826e-01 : 6.93e-07\n",
" 89 : 2.36e+01 : 3.55e-01 : 9.9999893113731242e-01 : 1.07e-06\n",
" 90 : 2.39e+01 : 3.22e-01 : 9.9999908937923931e-01 : 9.11e-07\n",
" 91 : 2.42e+01 : 3.00e-01 : 9.9999948312261555e-01 : 5.17e-07\n",
" 92 : 2.46e+01 : 3.14e-01 : 9.9999962343891935e-01 : 3.77e-07\n",
" 93 : 2.49e+01 : 3.49e-01 : 9.9999947771292474e-01 : 5.22e-07\n",
" 94 : 2.53e+01 : 3.65e-01 : 9.9999903846691773e-01 : 9.62e-07\n",
" 95 : 2.56e+01 : 3.04e-01 : 9.9999941461789033e-01 : 5.85e-07\n",
" 96 : 2.59e+01 : 3.04e-01 : 9.9999964361890326e-01 : 3.56e-07\n",
" 97 : 2.62e+01 : 3.42e-01 : 9.9999956280806190e-01 : 4.37e-07\n",
" 98 : 2.66e+01 : 3.70e-01 : 9.9999911836123156e-01 : 8.82e-07\n",
" 99 : 2.69e+01 : 3.08e-01 : 9.9999942637419115e-01 : 5.74e-07\n",
" 100 : 2.72e+01 : 3.08e-01 : 9.9999962677847165e-01 : 3.73e-07\n",
" 101 : 2.76e+01 : 3.43e-01 : 9.9999953595207036e-01 : 4.64e-07\n",
" 102 : 2.79e+01 : 3.67e-01 : 9.9999911916625339e-01 : 8.81e-07\n",
" 103 : 2.82e+01 : 3.08e-01 : 9.9999942036062439e-01 : 5.80e-07\n",
" 104 : 2.85e+01 : 3.08e-01 : 9.9999961856388120e-01 : 3.81e-07\n",
" 105 : 2.89e+01 : 3.42e-01 : 9.9999953272481978e-01 : 4.67e-07\n",
" 106 : 2.92e+01 : 3.65e-01 : 9.9999913473170166e-01 : 8.65e-07\n",
" 107 : 2.96e+01 : 3.09e-01 : 9.9999941835037198e-01 : 5.82e-07\n",
" 108 : 2.99e+01 : 3.09e-01 : 9.9999960900417784e-01 : 3.91e-07\n",
" 109 : 3.02e+01 : 3.41e-01 : 9.9999952481398846e-01 : 4.75e-07\n",
" 110 : 3.06e+01 : 3.64e-01 : 9.9999914583192606e-01 : 8.54e-07\n",
" 111 : 3.09e+01 : 3.10e-01 : 9.9999941566963535e-01 : 5.84e-07\n",
" 112 : 3.12e+01 : 3.10e-01 : 9.9999960026371204e-01 : 4.00e-07\n",
" 113 : 3.15e+01 : 3.41e-01 : 9.9999951840062606e-01 : 4.82e-07\n",
" 114 : 3.19e+01 : 3.62e-01 : 9.9999915670767547e-01 : 8.43e-07\n",
" 115 : 3.22e+01 : 3.45e-01 : 9.9999890421242110e-01 : 1.10e-06\n",
" 116 : 3.25e+01 : 3.11e-01 : 9.9999923763525689e-01 : 7.62e-07\n",
" 117 : 3.28e+01 : 3.00e-01 : 9.9999956557897374e-01 : 4.34e-07\n",
" 118 : 3.32e+01 : 3.25e-01 : 9.9999960936214261e-01 : 3.91e-07\n",
" 119 : 3.35e+01 : 3.59e-01 : 9.9999935128652762e-01 : 6.49e-07\n",
" 120 : 3.39e+01 : 3.60e-01 : 9.9999889546190690e-01 : 1.10e-06\n",
" 121 : 3.42e+01 : 3.25e-01 : 9.9999899783670898e-01 : 1.00e-06\n",
" 122 : 3.45e+01 : 2.98e-01 : 9.9999945697633119e-01 : 5.43e-07\n",
" 123 : 3.48e+01 : 3.08e-01 : 9.9999964335628166e-01 : 3.57e-07\n",
" 124 : 3.52e+01 : 3.46e-01 : 9.9999952622706911e-01 : 4.74e-07\n",
" 125 : 3.55e+01 : 3.69e-01 : 9.9999905886056595e-01 : 9.41e-07\n",
" 126 : 3.58e+01 : 3.04e-01 : 9.9999942282166476e-01 : 5.77e-07\n",
" 127 : 3.61e+01 : 3.04e-01 : 9.9999964603031355e-01 : 3.54e-07\n",
" 128 : 3.65e+01 : 3.43e-01 : 9.9999955870236679e-01 : 4.41e-07\n",
" 129 : 3.69e+01 : 3.71e-01 : 9.9999910499451694e-01 : 8.95e-07\n",
" 130 : 3.72e+01 : 3.07e-01 : 9.9999942649048967e-01 : 5.74e-07\n",
" 131 : 3.75e+01 : 3.07e-01 : 9.9999963250151580e-01 : 3.67e-07\n",
" 132 : 3.78e+01 : 3.43e-01 : 9.9999954173794070e-01 : 4.58e-07\n",
" 133 : 3.82e+01 : 3.68e-01 : 9.9999911352086013e-01 : 8.86e-07\n",
" 134 : 3.85e+01 : 3.08e-01 : 9.9999942207296633e-01 : 5.78e-07\n",
" 135 : 3.88e+01 : 3.08e-01 : 9.9999962322897262e-01 : 3.77e-07\n",
" 136 : 3.91e+01 : 3.42e-01 : 9.9999953599883418e-01 : 4.64e-07\n",
" 137 : 3.95e+01 : 3.66e-01 : 9.9999912820749648e-01 : 8.72e-07\n",
" 138 : 3.98e+01 : 3.09e-01 : 9.9999941961677052e-01 : 5.80e-07\n",
" 139 : 4.01e+01 : 3.09e-01 : 9.9999961361827316e-01 : 3.86e-07\n",
" 140 : 4.05e+01 : 3.42e-01 : 9.9999952836548056e-01 : 4.72e-07\n",
" 141 : 4.08e+01 : 3.64e-01 : 9.9999914014091940e-01 : 8.60e-07\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" 142 : 4.11e+01 : 3.10e-01 : 9.9999941693883510e-01 : 5.83e-07\n",
" 143 : 4.14e+01 : 3.10e-01 : 9.9999960463251514e-01 : 3.95e-07\n",
" 144 : 4.18e+01 : 3.41e-01 : 9.9999952165610273e-01 : 4.78e-07\n",
" 145 : 4.22e+01 : 3.63e-01 : 9.9999915140461471e-01 : 8.49e-07\n",
" 146 : 4.25e+01 : 3.11e-01 : 9.9999941448747354e-01 : 5.86e-07\n",
" 147 : 4.28e+01 : 3.11e-01 : 9.9999959600897592e-01 : 4.04e-07\n",
" 148 : 4.31e+01 : 3.41e-01 : 9.9999951520446739e-01 : 4.85e-07\n",
" 149 : 4.35e+01 : 3.61e-01 : 9.9999916172552206e-01 : 8.38e-07\n",
" 150 : 4.38e+01 : 3.45e-01 : 9.9999891720742473e-01 : 1.08e-06\n",
" 151 : 4.41e+01 : 3.12e-01 : 9.9999924067811397e-01 : 7.59e-07\n",
" 152 : 4.44e+01 : 3.01e-01 : 9.9999956202100138e-01 : 4.38e-07\n",
" 153 : 4.48e+01 : 3.25e-01 : 9.9999960508702357e-01 : 3.95e-07\n",
" 154 : 4.51e+01 : 3.58e-01 : 9.9999935116023642e-01 : 6.49e-07\n",
" 155 : 4.55e+01 : 3.60e-01 : 9.9999890813343284e-01 : 1.09e-06\n",
" 156 : 4.58e+01 : 3.25e-01 : 9.9999900756889537e-01 : 9.92e-07\n",
" 157 : 4.61e+01 : 2.98e-01 : 9.9999945545984126e-01 : 5.45e-07\n",
" 158 : 4.64e+01 : 3.09e-01 : 9.9999963905838241e-01 : 3.61e-07\n",
" 159 : 4.68e+01 : 3.46e-01 : 9.9999952293759720e-01 : 4.77e-07\n",
" 160 : 4.71e+01 : 3.68e-01 : 9.9999906568072106e-01 : 9.34e-07\n",
" 161 : 4.74e+01 : 3.05e-01 : 9.9999942167245681e-01 : 5.78e-07\n",
" 162 : 4.77e+01 : 3.05e-01 : 9.9999964202521041e-01 : 3.58e-07\n",
" 163 : 4.81e+01 : 3.43e-01 : 9.9999955552367414e-01 : 4.44e-07\n",
" 164 : 4.84e+01 : 3.70e-01 : 9.9999911060056901e-01 : 8.89e-07\n",
" 165 : 4.88e+01 : 3.07e-01 : 9.9999942534481068e-01 : 5.75e-07\n",
" 166 : 4.91e+01 : 3.07e-01 : 9.9999962870609627e-01 : 3.71e-07\n",
" 167 : 4.94e+01 : 3.43e-01 : 9.9999953881002768e-01 : 4.61e-07\n",
" 168 : 4.98e+01 : 3.67e-01 : 9.9999911871232416e-01 : 8.81e-07\n",
" 169 : 5.01e+01 : 3.08e-01 : 9.9999942101613537e-01 : 5.79e-07\n",
" 170 : 5.04e+01 : 3.08e-01 : 9.9999961962214523e-01 : 3.80e-07\n",
" 171 : 5.07e+01 : 3.42e-01 : 9.9999953323688884e-01 : 4.67e-07\n",
" 172 : 5.11e+01 : 3.65e-01 : 9.9999913291306763e-01 : 8.67e-07\n",
" 173 : 5.14e+01 : 3.09e-01 : 9.9999941861799835e-01 : 5.81e-07\n",
" 174 : 5.17e+01 : 3.09e-01 : 9.9999961018322492e-01 : 3.90e-07\n",
" 175 : 5.21e+01 : 3.41e-01 : 9.9999952576749485e-01 : 4.74e-07\n",
" 176 : 5.24e+01 : 3.64e-01 : 9.9999914445875404e-01 : 8.56e-07\n",
" 177 : 5.27e+01 : 3.10e-01 : 9.9999941599720532e-01 : 5.84e-07\n",
" 178 : 5.30e+01 : 3.10e-01 : 9.9999960135263399e-01 : 3.99e-07\n",
" 179 : 5.34e+01 : 3.41e-01 : 9.9999951920016095e-01 : 4.81e-07\n",
" 180 : 5.37e+01 : 3.62e-01 : 9.9999915537954931e-01 : 8.45e-07\n",
" 181 : 5.41e+01 : 3.45e-01 : 9.9999890085620846e-01 : 1.10e-06\n",
" 182 : 5.44e+01 : 3.11e-01 : 9.9999923688555781e-01 : 7.63e-07\n",
" 183 : 5.47e+01 : 3.00e-01 : 9.9999956649957022e-01 : 4.34e-07\n",
" 184 : 5.50e+01 : 3.25e-01 : 9.9999961044319396e-01 : 3.90e-07\n",
" 185 : 5.54e+01 : 3.59e-01 : 9.9999935128701556e-01 : 6.49e-07\n",
" 186 : 5.57e+01 : 3.61e-01 : 9.9999889216712323e-01 : 1.11e-06\n",
" 187 : 5.61e+01 : 3.25e-01 : 9.9999899536050973e-01 : 1.00e-06\n",
" 188 : 5.64e+01 : 2.98e-01 : 9.9999945738735840e-01 : 5.43e-07\n",
" 189 : 5.67e+01 : 3.08e-01 : 9.9999964444946321e-01 : 3.56e-07\n",
" 190 : 5.70e+01 : 3.46e-01 : 9.9999952704483852e-01 : 4.73e-07\n",
" 191 : 5.74e+01 : 3.70e-01 : 9.9999905706846193e-01 : 9.43e-07\n",
" 192 : 5.77e+01 : 3.04e-01 : 9.9999942311229106e-01 : 5.77e-07\n",
" 193 : 5.80e+01 : 3.04e-01 : 9.9999964705875743e-01 : 3.53e-07\n",
" 194 : 5.83e+01 : 3.43e-01 : 9.9999955952555675e-01 : 4.40e-07\n",
" 195 : 5.87e+01 : 3.71e-01 : 9.9999910354638777e-01 : 8.96e-07\n",
" 196 : 5.90e+01 : 3.07e-01 : 9.9999942678583198e-01 : 5.73e-07\n",
" 197 : 5.93e+01 : 3.07e-01 : 9.9999963347296739e-01 : 3.67e-07\n",
" 198 : 5.97e+01 : 3.43e-01 : 9.9999954248858458e-01 : 4.58e-07\n",
" 199 : 6.00e+01 : 3.68e-01 : 9.9999911217688986e-01 : 8.88e-07\n",
" 200 : 6.03e+01 : 3.08e-01 : 9.9999942234397543e-01 : 5.78e-07\n",
" 201 : 6.06e+01 : 3.08e-01 : 9.9999962415206511e-01 : 3.76e-07\n",
" 202 : 6.10e+01 : 3.42e-01 : 9.9999953670784880e-01 : 4.63e-07\n",
" 203 : 6.14e+01 : 3.66e-01 : 9.9999912699211180e-01 : 8.73e-07\n",
" 204 : 6.17e+01 : 3.09e-01 : 9.9999941987296970e-01 : 5.80e-07\n",
" 205 : 6.20e+01 : 3.09e-01 : 9.9999961449675778e-01 : 3.86e-07\n",
" 206 : 6.23e+01 : 3.42e-01 : 9.9999952903153067e-01 : 4.71e-07\n",
" 207 : 6.27e+01 : 3.64e-01 : 9.9999913902682047e-01 : 8.61e-07\n",
" 208 : 6.30e+01 : 3.10e-01 : 9.9999941718016094e-01 : 5.83e-07\n",
" 209 : 6.33e+01 : 3.10e-01 : 9.9999960547091049e-01 : 3.95e-07\n",
" 210 : 6.36e+01 : 3.41e-01 : 9.9999952228537503e-01 : 4.78e-07\n",
" 211 : 6.40e+01 : 3.63e-01 : 9.9999915038014564e-01 : 8.50e-07\n",
" 212 : 6.43e+01 : 3.11e-01 : 9.9999941471576614e-01 : 5.85e-07\n",
" 213 : 6.46e+01 : 3.11e-01 : 9.9999959681069994e-01 : 4.03e-07\n",
" 214 : 6.50e+01 : 3.41e-01 : 9.9999951580028834e-01 : 4.84e-07\n",
" 215 : 6.53e+01 : 3.61e-01 : 9.9999916077889850e-01 : 8.39e-07\n",
" 216 : 6.57e+01 : 3.45e-01 : 9.9999891478617053e-01 : 1.09e-06\n",
" 217 : 6.60e+01 : 3.11e-01 : 9.9999924011841390e-01 : 7.60e-07\n",
" 218 : 6.63e+01 : 3.01e-01 : 9.9999956269483181e-01 : 4.37e-07\n",
" 219 : 6.66e+01 : 3.25e-01 : 9.9999960589031367e-01 : 3.94e-07\n",
" 220 : 6.70e+01 : 3.58e-01 : 9.9999935117444894e-01 : 6.49e-07\n",
" 221 : 6.73e+01 : 3.60e-01 : 9.9999890576137107e-01 : 1.09e-06\n",
" 222 : 6.77e+01 : 3.25e-01 : 9.9999900576076706e-01 : 9.94e-07\n",
" 223 : 6.80e+01 : 2.98e-01 : 9.9999945575384608e-01 : 5.44e-07\n",
" 224 : 6.83e+01 : 3.08e-01 : 9.9999963987041862e-01 : 3.60e-07\n",
" 225 : 6.86e+01 : 3.46e-01 : 9.9999952354994859e-01 : 4.76e-07\n",
" 226 : 6.90e+01 : 3.69e-01 : 9.9999906438984965e-01 : 9.36e-07\n",
" 227 : 6.93e+01 : 3.05e-01 : 9.9999942188864188e-01 : 5.78e-07\n",
" 228 : 6.96e+01 : 3.05e-01 : 9.9999964278632258e-01 : 3.57e-07\n",
" 229 : 6.99e+01 : 3.43e-01 : 9.9999955612794800e-01 : 4.44e-07\n",
" 230 : 7.03e+01 : 3.70e-01 : 9.9999910954510063e-01 : 8.90e-07\n",
" 231 : 7.06e+01 : 3.07e-01 : 9.9999942556231969e-01 : 5.74e-07\n",
" 232 : 7.09e+01 : 3.07e-01 : 9.9999962942688203e-01 : 3.71e-07\n",
" 233 : 7.13e+01 : 3.43e-01 : 9.9999953936468555e-01 : 4.61e-07\n",
" 234 : 7.16e+01 : 3.67e-01 : 9.9999911773236638e-01 : 8.82e-07\n",
" 235 : 7.19e+01 : 3.08e-01 : 9.9999942121652141e-01 : 5.79e-07\n",
" 236 : 7.22e+01 : 3.08e-01 : 9.9999962030760037e-01 : 3.80e-07\n",
" 237 : 7.26e+01 : 3.42e-01 : 9.9999953376092610e-01 : 4.66e-07\n",
" 238 : 7.29e+01 : 3.66e-01 : 9.9999913202433444e-01 : 8.68e-07\n",
" 239 : 7.33e+01 : 3.09e-01 : 9.9999941880754306e-01 : 5.81e-07\n",
" 240 : 7.36e+01 : 3.09e-01 : 9.9999961083624178e-01 : 3.89e-07\n",
" 241 : 7.39e+01 : 3.41e-01 : 9.9999952626056499e-01 : 4.74e-07\n",
" 242 : 7.43e+01 : 3.64e-01 : 9.9999914364254705e-01 : 8.56e-07\n",
" 243 : 7.46e+01 : 3.10e-01 : 9.9999941617596422e-01 : 5.84e-07\n",
" 244 : 7.49e+01 : 3.10e-01 : 9.9999960197636661e-01 : 3.98e-07\n",
" 245 : 7.52e+01 : 3.41e-01 : 9.9999951966650213e-01 : 4.80e-07\n",
" 246 : 7.56e+01 : 3.62e-01 : 9.9999915462765532e-01 : 8.45e-07\n",
" 247 : 7.59e+01 : 3.45e-01 : 9.9999889891782989e-01 : 1.10e-06\n",
" 248 : 7.62e+01 : 3.11e-01 : 9.9999923644016553e-01 : 7.64e-07\n",
" 249 : 7.65e+01 : 3.00e-01 : 9.9999956702241355e-01 : 4.33e-07\n",
" 250 : 7.69e+01 : 3.25e-01 : 9.9999961106656898e-01 : 3.89e-07\n",
" 251 : 7.72e+01 : 3.59e-01 : 9.9999935130057604e-01 : 6.49e-07\n",
" 252 : 7.76e+01 : 3.61e-01 : 9.9999889027711686e-01 : 1.11e-06\n",
" 253 : 7.79e+01 : 3.25e-01 : 9.9999899391941260e-01 : 1.01e-06\n",
" 254 : 7.82e+01 : 2.97e-01 : 9.9999945761217401e-01 : 5.42e-07\n",
" 255 : 7.85e+01 : 3.08e-01 : 9.9999964507518202e-01 : 3.55e-07\n",
" 256 : 7.89e+01 : 3.46e-01 : 9.9999952752346555e-01 : 4.72e-07\n",
" 257 : 7.92e+01 : 3.70e-01 : 9.9999905605503181e-01 : 9.44e-07\n",
" 258 : 7.95e+01 : 3.04e-01 : 9.9999942327967739e-01 : 5.77e-07\n",
" 259 : 7.98e+01 : 3.04e-01 : 9.9999964764224414e-01 : 3.52e-07\n",
" 260 : 8.02e+01 : 3.43e-01 : 9.9999955999117707e-01 : 4.40e-07\n",
" 261 : 8.06e+01 : 3.71e-01 : 9.9999910271871828e-01 : 8.97e-07\n",
" 262 : 8.09e+01 : 3.07e-01 : 9.9999942695334310e-01 : 5.73e-07\n",
" 263 : 8.12e+01 : 3.07e-01 : 9.9999963402505132e-01 : 3.66e-07\n",
" 264 : 8.15e+01 : 3.43e-01 : 9.9999954291601634e-01 : 4.57e-07\n",
" 265 : 8.19e+01 : 3.68e-01 : 9.9999911141127928e-01 : 8.89e-07\n",
" 266 : 8.22e+01 : 3.07e-01 : 9.9999942249813900e-01 : 5.78e-07\n",
" 267 : 8.25e+01 : 3.07e-01 : 9.9999962467630787e-01 : 3.75e-07\n",
" 268 : 8.28e+01 : 3.42e-01 : 9.9999953711077494e-01 : 4.63e-07\n",
" 269 : 8.32e+01 : 3.66e-01 : 9.9999912629963417e-01 : 8.74e-07\n",
" 270 : 8.35e+01 : 3.09e-01 : 9.9999942001856534e-01 : 5.80e-07\n",
" 271 : 8.38e+01 : 3.09e-01 : 9.9999961499562351e-01 : 3.85e-07\n",
" 272 : 8.42e+01 : 3.42e-01 : 9.9999952941008663e-01 : 4.71e-07\n",
" 273 : 8.45e+01 : 3.65e-01 : 9.9999913839239463e-01 : 8.62e-07\n",
" 274 : 8.48e+01 : 3.10e-01 : 9.9999941731729769e-01 : 5.83e-07\n",
" 275 : 8.52e+01 : 3.10e-01 : 9.9999960594691883e-01 : 3.94e-07\n",
" 276 : 8.55e+01 : 3.41e-01 : 9.9999952264291714e-01 : 4.77e-07\n",
" 277 : 8.59e+01 : 3.63e-01 : 9.9999914979694582e-01 : 8.50e-07\n",
" 278 : 8.62e+01 : 3.11e-01 : 9.9999941484546684e-01 : 5.85e-07\n",
" 279 : 8.65e+01 : 3.11e-01 : 9.9999959726582233e-01 : 4.03e-07\n",
" 280 : 8.68e+01 : 3.41e-01 : 9.9999951613876403e-01 : 4.84e-07\n",
" 281 : 8.72e+01 : 3.61e-01 : 9.9999916024018709e-01 : 8.40e-07\n",
" 282 : 8.75e+01 : 3.45e-01 : 9.9999891340699587e-01 : 1.09e-06\n",
" 283 : 8.78e+01 : 3.11e-01 : 9.9999923979983552e-01 : 7.60e-07\n",
" 284 : 8.81e+01 : 3.01e-01 : 9.9999956307723359e-01 : 4.37e-07\n",
" 285 : 8.85e+01 : 3.25e-01 : 9.9999960634618879e-01 : 3.94e-07\n",
" 286 : 8.88e+01 : 3.59e-01 : 9.9999935118272609e-01 : 6.49e-07\n",
" 287 : 8.92e+01 : 3.60e-01 : 9.9999890441095918e-01 : 1.10e-06\n",
" 288 : 8.95e+01 : 3.25e-01 : 9.9999900473138170e-01 : 9.95e-07\n",
" 289 : 8.98e+01 : 2.98e-01 : 9.9999945592041528e-01 : 5.44e-07\n",
" 290 : 9.01e+01 : 3.08e-01 : 9.9999964033087207e-01 : 3.60e-07\n",
" 291 : 9.05e+01 : 3.46e-01 : 9.9999952389776126e-01 : 4.76e-07\n",
" 292 : 9.08e+01 : 3.69e-01 : 9.9999906365622926e-01 : 9.36e-07\n",
" 293 : 9.11e+01 : 3.05e-01 : 9.9999942201129732e-01 : 5.78e-07\n",
" 294 : 9.14e+01 : 3.05e-01 : 9.9999964321763968e-01 : 3.57e-07\n",
" 295 : 9.18e+01 : 3.43e-01 : 9.9999955647059524e-01 : 4.44e-07\n",
" 296 : 9.21e+01 : 3.70e-01 : 9.9999910894535204e-01 : 8.91e-07\n",
" 297 : 9.25e+01 : 3.07e-01 : 9.9999942568564926e-01 : 5.74e-07\n",
" 298 : 9.28e+01 : 3.07e-01 : 9.9999962983530322e-01 : 3.70e-07\n",
" 299 : 9.31e+01 : 3.43e-01 : 9.9999953967920163e-01 : 4.60e-07\n",
" 300 : 9.35e+01 : 3.67e-01 : 9.9999911717576606e-01 : 8.83e-07\n",
" 301 : 9.38e+01 : 3.08e-01 : 9.9999942133012731e-01 : 5.79e-07\n",
" 302 : 9.41e+01 : 3.08e-01 : 9.9999962069593407e-01 : 3.79e-07\n",
" 303 : 9.44e+01 : 3.42e-01 : 9.9999953405800079e-01 : 4.66e-07\n",
" 304 : 9.48e+01 : 3.66e-01 : 9.9999913151971154e-01 : 8.68e-07\n",
" 305 : 9.51e+01 : 3.09e-01 : 9.9999941891498323e-01 : 5.81e-07\n",
" 306 : 9.54e+01 : 3.09e-01 : 9.9999961120614789e-01 : 3.89e-07\n",
" 307 : 9.58e+01 : 3.41e-01 : 9.9999952654003588e-01 : 4.73e-07\n",
" 308 : 9.61e+01 : 3.64e-01 : 9.9999914317923666e-01 : 8.57e-07\n",
" 309 : 9.64e+01 : 3.10e-01 : 9.9999941627727518e-01 : 5.84e-07\n",
" 310 : 9.67e+01 : 3.10e-01 : 9.9999960232964236e-01 : 3.98e-07\n",
" 311 : 9.71e+01 : 3.41e-01 : 9.9999951993077896e-01 : 4.80e-07\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" 312 : 9.74e+01 : 3.62e-01 : 9.9999915420095764e-01 : 8.46e-07\n",
" 313 : 9.78e+01 : 3.45e-01 : 9.9999889781701456e-01 : 1.10e-06\n",
" 314 : 9.81e+01 : 3.11e-01 : 9.9999923618737285e-01 : 7.64e-07\n",
" 315 : 9.84e+01 : 3.00e-01 : 9.9999956731847128e-01 : 4.33e-07\n",
" 316 : 9.87e+01 : 3.25e-01 : 9.9999961141955540e-01 : 3.89e-07\n",
" 317 : 9.91e+01 : 3.59e-01 : 9.9999935130838857e-01 : 6.49e-07\n",
" 318 : 9.94e+01 : 3.61e-01 : 9.9999888920425595e-01 : 1.11e-06\n",
" 319 : 9.98e+01 : 3.25e-01 : 9.9999899310135654e-01 : 1.01e-06\n",
" 320 : 1.00e+02 : 2.30e-01 : 9.9999980475946604e-01 : 1.95e-07\n"
]
}
],
"source": [
"println(\"step : h : t : y : error\")\n",
"\n",
"stepcnt = 0\n",
"while sol.status != \"finished\"\n",
" sol.step()\n",
" global stepcnt = stepcnt + 1\n",
" scnt = @sprintf(\"%4d\", stepcnt)\n",
" step = @sprintf(\"%.2e\", sol.step_size)\n",
" tval = @sprintf(\"%.2e\", sol.t)\n",
" soly = @sprintf(\"%.16e\", sol.y[1])\n",
" exact = 1 - exp(-10*sol.t)/2\n",
" err = abs(sol.y[1] - exact)\n",
" serr = @sprintf(\"%.2e\", err)\n",
" println(\"$scnt : $tval : $step : $soly : $serr\")\n",
"end"
]
},
{
"cell_type": "markdown",
"id": "6817325f",
"metadata": {},
"source": [
"We see that ``RK45`` of the ``integrate`` module of ``SciPy`` requires 320 steps."
]
},
{
"cell_type": "markdown",
"id": "329bc179",
"metadata": {},
"source": [
"Let us run the backward Euler method on the same problem."
]
},
{
"cell_type": "code",
"execution_count": 24,
"id": "4cb21ada",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"backwardEuler (generic function with 1 method)"
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function backwardEuler(n::Int)\n",
" h = 100.0/n\n",
" y0 = 0.5\n",
" println(\"step : x : approximation : error\")\n",
" for i=1:n\n",
" x = i*h\n",
" y1 = (y0 + 10*h)/(1 + 10*h)\n",
" y0 = y1\n",
" exact = 1 - exp(-10*x)/2\n",
" err = abs(y1 - exact)\n",
" scnt = @sprintf(\"%4d\", i)\n",
" strx = @sprintf(\"%.2e\", x)\n",
" ssol = @sprintf(\"%.16e\", y1)\n",
" serr = @sprintf(\"%.2e\", err)\n",
" println(\"$scnt : $strx : $ssol : $serr\")\n",
" end\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 25,
"id": "eef68ed6",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"step : x : approximation : error\n",
" 1 : 1.00e+01 : 9.9504950495049505e-01 : 4.95e-03\n",
" 2 : 2.00e+01 : 9.9995098519752967e-01 : 4.90e-05\n",
" 3 : 3.00e+01 : 9.9999951470492598e-01 : 4.85e-07\n",
" 4 : 4.00e+01 : 9.9999999519509830e-01 : 4.80e-09\n",
" 5 : 5.00e+01 : 9.9999999995242672e-01 : 4.76e-11\n",
" 6 : 6.00e+01 : 9.9999999999952893e-01 : 4.71e-13\n",
" 7 : 7.00e+01 : 9.9999999999999534e-01 : 4.66e-15\n",
" 8 : 8.00e+01 : 1.0000000000000000e+00 : 0.00e+00\n",
" 9 : 9.00e+01 : 1.0000000000000000e+00 : 0.00e+00\n",
" 10 : 1.00e+02 : 1.0000000000000000e+00 : 0.00e+00\n"
]
}
],
"source": [
"backwardEuler(10)"
]
},
{
"cell_type": "markdown",
"id": "3523f86b",
"metadata": {},
"source": [
"We see that with only 10 steps, we achieve rapid convergence."
]
},
{
"cell_type": "markdown",
"id": "8f8115c5",
"metadata": {},
"source": [
"## 4. Using DifferentialEquations"
]
},
{
"cell_type": "markdown",
"id": "e5328f03",
"metadata": {},
"source": [
"Consider the solution of\n",
"\n",
"$$\n",
" \\frac{d}{d~\\!t} y = y, \\quad y(0) = y_0, \\quad t \\in [0,1]\n",
"$$\n",
"\n",
"with a fourth order Runge-Kutta method."
]
},
{
"cell_type": "code",
"execution_count": 26,
"id": "8d719811",
"metadata": {},
"outputs": [],
"source": [
"using DifferentialEquations"
]
},
{
"cell_type": "code",
"execution_count": 27,
"id": "f1c7c868",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"testfun (generic function with 1 method)"
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"testfun(y, p, t) = y"
]
},
{
"cell_type": "code",
"execution_count": 28,
"id": "7f221935",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1"
]
},
"execution_count": 28,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"y0 = 1"
]
},
{
"cell_type": "code",
"execution_count": 29,
"id": "0b2fce4a",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(0.0, 1.0)"
]
},
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"timespan = (0.0, 1.0)"
]
},
{
"cell_type": "code",
"execution_count": 30,
"id": "108a76bf",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"\u001b[36mODEProblem\u001b[0m with uType \u001b[36mInt64\u001b[0m and tType \u001b[36mFloat64\u001b[0m. In-place: \u001b[36mfalse\u001b[0m\n",
"timespan: (0.0, 1.0)\n",
"u0: 1"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ode = ODEProblem(testfun, y0, timespan)"
]
},
{
"cell_type": "code",
"execution_count": 31,
"id": "d1180217",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"retcode: Success\n",
"Interpolation: 3rd order Hermite\n",
"t: 53-element Vector{Float64}:\n",
" 0.0\n",
" 0.0037606030930863936\n",
" 0.009904232271738917\n",
" 0.01741266492992194\n",
" 0.02711332433300147\n",
" 0.0384578161191669\n",
" 0.05162397134853269\n",
" 0.06627222802749277\n",
" 0.08231516797356161\n",
" 0.09950872249018991\n",
" 0.11770152088033474\n",
" 0.1367114189918031\n",
" 0.1563998303815034\n",
" ⋮\n",
" 0.7840950043865613\n",
" 0.8055204080644183\n",
" 0.8269095529529708\n",
" 0.8482629408911011\n",
" 0.8695810911470361\n",
" 0.8908645403256985\n",
" 0.9121138352555783\n",
" 0.9333295329494489\n",
" 0.9545121851605235\n",
" 0.9756623538268859\n",
" 0.9967805900913862\n",
" 1.0\n",
"u: 53-element Vector{Float64}:\n",
" 1.0\n",
" 1.0037676830330582\n",
" 1.0099534415058256\n",
" 1.0175651491460624\n",
" 1.0274842351306568\n",
" 1.039206889651116\n",
" 1.0529797175103524\n",
" 1.0685175579917117\n",
" 1.0857979646309999\n",
" 1.1046281057699228\n",
" 1.1249082995595794\n",
" 1.146497243366542\n",
" 1.1692936286466074\n",
" ⋮\n",
" 2.19042371692509\n",
" 2.2378607938210835\n",
" 2.286242297480796\n",
" 2.3355862716854854\n",
" 2.3859111615701503\n",
" 2.4372358259508022\n",
" 2.4895795326139774\n",
" 2.5429619698018335\n",
" 2.5974032180629654\n",
" 2.6529237993061145\n",
" 2.709544632569373\n",
" 2.718281824161426"
]
},
"execution_count": 31,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sol = solve(ode, RK4(), reltol=1e-8, abstol=1e-8)"
]
},
{
"cell_type": "code",
"execution_count": 32,
"id": "a33d7470",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 32,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(sol.t, sol.u, label=\"sol\", legend=:topleft)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.8.0",
"language": "julia",
"name": "julia-1.8"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.8.0"
},
"latex_envs": {
"LaTeX_envs_menu_present": true,
"autocomplete": true,
"bibliofile": "biblio.bib",
"cite_by": "apalike",
"current_citInitial": 1,
"eqLabelWithNumbers": true,
"eqNumInitial": 1,
"hotkeys": {
"equation": "Ctrl-E",
"itemize": "Ctrl-I"
},
"labels_anchors": false,
"latex_user_defs": false,
"report_style_numbering": false,
"user_envs_cfg": false
}
},
"nbformat": 4,
"nbformat_minor": 5
}