{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Numerisk integrasjon\n",
"** Anne Kværnø **"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Problemstilling\n",
"\n",
"Vi skal altså finne en numerisk tilnærmelse til integralet\n",
"$$ I(a,b) = \\int_a^b f(x) dx $$\n",
"for en gitt funksjon $f(x)$. \n",
"\n",
"Teknikken vi skal diskutere kalles _numeriske kvadratur_, disse er på formen\n",
"\n",
"$$ Q(a,b) = \\sum_{i=0}^n w_i f(x_i) $$\n",
"\n",
"der punktene $x_i$ kalles noder og $w_i$ vekter. Og vi ønsker altså at \n",
"\n",
"$$ I(a,b) \\approx Q(a,b). $$\n",
"\n",
"[Trapesregelen, midtpunkt-regelen og Simpsons regel](https://wiki.math.ntnu.no/tma4100/tema/numerics?numerisk_integrasjon)\n",
"er alle kjente eksempler på slike formler. \n",
"\n",
"I dette notatet skal vi se mer generelt på hvordan slike kvadraturformler kan konstrueres, hvordan man kan finne et estimat for feilen i kvadratueret, og hvordan dette kan settes sammen til en algoritme som gir et tilnærmelse til løsningen med en nøyaktighet som er mindre eller omtrent lik en oppgitt toleranse. \n",
"\n",
"\n",
"# Kvadraturformler basert på polynominterpolasjon. \n",
"Velg et sett med noder $x_i$, $i=0,\\dotsc,n$ på intervallet $[a,b]$. Interpolasjonspolynomet til $f(x)$ gjennom disse nodene er gitt ved\n",
"\n",
"$$ p_n(x) = \\sum_{i=0}^n f(x_i) \\ell_i(x) $$ \n",
"\n",
"der _kardinalfunksjonene_ $\\ell_i(x)$ er gitt ved \n",
"\n",
"$$\\ell_i(x) = \\prod_{j=0,j\\not=i}^n \\frac{x-x_j}{x_i-x_j},$$\n",
"\n",
"se notatet om polynominterpolasjon. Da kan vi bruke integralet av $p_n(x)$ som en tilnærmelse til integralet av $f(x)$, dvs\n",
"\n",
"$$\n",
"I(a,b) = \\int_{a}^b f(x)dx \\approx \\int_a^b p_n(x)dx \n",
" = \\sum_{i=0}^n f(x_i) \\int_a^b \\ell_i(x) dx = \\sum_{i=0}^n w_i f(x_i) = Q(a,b) \n",
"$$\n",
"\n",
"##### Definisjon: Presisjonsgrad\n",
"Et numerisk kvadratur har presisjonsgrad $d$ dersom \n",
"$Q(a,b) = I(a,b)$ for alle polynomer av grad $d$ eller mindre. \n",
"\n",
"\n",
"\n",
"##### Eksempler\n",
"1. Trapesmetoden er basert på nodene $(a,b)$, og blir\n",
" $$ T(a,b) = \\frac{b-a}{2}(f(a)+f(b)). $$\n",
" Presisjonsgraden finner vi ved: \n",
" \n",
" \\begin{align}\n",
" f &= 1, & I(a,b) &= b-a, & T(a,b) &= b-a, \\\\\n",
" f &= x, & I(a,b) &= \\frac{b^2-a^2}{2}, & T(a,b) &= \\frac{b^2-a^2}{2}, \\\\\n",
" f &= x^2, & I(a,b) &= \\frac{b^3-a^3}{3}, & T(a,b) &= \\frac{(b-a)(b^2+a^2)}{2}, \n",
" \\end{align}\n",
" \n",
" så trapesmetoden har presisjonsgrad 1. \n",
"\n",
"2. På intervallet $[-1,1]$, velg noder $\\pm \\sqrt(3)/3$. Dette gir kvadraturformelen \n",
" $$ Q(-1, 1) = f\\left(-\\frac{\\sqrt{3}}{3}\\right) + f\\left(\\frac{\\sqrt{3}}{3}\\right). $$\n",
" Denne har har også presisjonsgrad 3 (vis det selv). \n",
" \n",
"Jo høyere presisjonsgrad, jo bedre kvadratur. \n",
"Alle kvadraturformler basert på polynominterpolasjon i $n+1$ noder vil automatisk ha presisjonsgrad minst $n$, men bedre resultat er som vi har sett mulig. \n",
"\n",
"I praksis vil vi utvikle slike formler på et standardinterval $[-1,1]$ for deretter å flytte disse over til intervallet $[a,b]$. Eller vi bruker _sammensatte formler_: Partisjoner intervallet $[a,b]$ i $N$ passende delintervaller\n",
"\n",
"$$ a = X_0 < X_1 < \\dotsm < X_N = b $$\n",
"\n",
"og anvend en kvadraturformel på hvert delintervall:\n",
"\n",
"$$\n",
"\\int_a^b f(x) = \\sum_{k=0}^{N-1} \\int_{X_k}^{X_{k+1}} f(x) dx \n",
"\\approx \\sum_{k=0}^{N-1} Q(X_{k}, X_{k+1}).\n",
"$$\n",
"\n",
"Oppsummert: \n",
"1. Velg noder på standardintervallet $[-1, 1]$. \n",
"2. Finn kvadraturformelen $Q(-1,1)$ ved å integrere interpolasjonspolynomet. \n",
"3. Flytt dette til et tilfeldig intervall $[a,b] \\; \\Rightarrow \\; Q(a,b)$. \n",
"4. Forutsatt at feilen $E(a,b) = I(a,b)-Q(a,b) \\approx C (b-a)^p$ for en gitt $p$ er det også mulig å finne et feilestimat. Det kan brukes til å konstruere en algoritme som automatisk gir en passende partisjonering av intervallet. \n",
"\n",
"La oss gå gjennom hele denne prosessen med et velkjent eksempel:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Simpsons metode\n",
"\n",
"##### Standardintervall [-1,1]:\n",
"\n",
"Vi starter med å finne en tilnærmelse til $ \\int_{-1}^1 f(t) dt$. \n",
"Velg nodene $t_0=-1$, $t_1=0$ og $t_2=1$. Dette gir oss kardinalfunksjonene\n",
"\n",
"$$ \n",
"\\ell_0 = \\frac{1}{2}(t^2-t), \\qquad \n",
"\\ell_1(t) = 1-t^2, \\qquad \n",
"\\ell_2(t) = \\frac{1}{2}(t^2+t). \n",
"$$\n",
"\n",
"som igjen gir vektene \n",
"\n",
"$$ \n",
"w_0 = \\int_{-1}^1 \\ell_0(t)dt = \\frac{1}{3}, \\qquad\n",
"w_1 = \\int_{-1}^1 \\ell_1(t)dt = \\frac{4}{3}, \\qquad\n",
"w_2 = \\int_{-1}^1 \\ell_2(t)dt = \\frac{1}{3}\n",
"$$\n",
"\n",
"slik at\n",
"\n",
"$$ \n",
"\\int_{-1}^1 f(t) dt \\approx \\int_{-1}^1 p_2(t) dt = \\sum_{i=0}^2 w_i f(t_i) = \n",
"\\frac{1}{3} \\left[\\; f(-1) + 4 f(0) + f(1) \\; \\right].\n",
"$$\n",
"\n",
"Simpsons metode har presisjonsgrad 3 (sjekk det selv). \n",
"\n",
"_Eksempel_: \n",
"$$ \n",
"\\int_{-1}^1 \\cos \\left( \\frac{\\pi t}{2}\\right)dt = \\frac{4}{\\pi} \n",
"\\approx \\frac{1}{3}\\left[\\cos(-\\pi/2) + 4 \\cos(0) \n",
"+ \\cos(\\pi/2) \\right]= \\frac{4}{3}. \n",
"$$\n",
"\n",
"###### Over til intervallet [a,b]:\n",
"\n",
"Resultatet fra standardintervallet flyttes til et intervall $[a,b]$ ved transformasjonen\n",
"\n",
"$$ x = \\frac{b-a}{2}t + \\frac{b+a}{2}. $$\n",
"\n",
"Husk at da blir $dx = (b-a)/2\\; dt$. \n",
"\n",
"Det betyr at \n",
"\n",
"$$\n",
"\\int_a^b f(x)dx = \\frac{b-a}{2} \\int_{-1}^1 f\\left(\\frac{b-a}{2}t + \\frac{b+a}{2}\\right) dt\n",
"\\approx \\frac{b-a}{6} \\left[\\;f(a)+4f\\left(\\frac{b+a}{2}\\right)+f(b)\\;\\right]. \n",
"$$\n",
"\n",
"Så Simpsons formel over et interval $[a,b]$ er\n",
"\n",
"$$ S(a,b) = \\frac{b-a}{6}\\left(\\; f(a)+4f(c)+f(b)\\; \\right), \\qquad c=\\frac{b+a}{2}. $$\n",
"\n",
"(Det er vanlig praksis å bruke $S$ i stedet for $Q$ i dette tilfellet.)\n",
"\n",
"###### Sammensatt Simpsons metode\n",
"Del intervallet $[a,b]$ inn i $2m$ like store intervaller av lengde \n",
"$ h = (b-a)/(2m)$. La $x_j = a+jh$, $i=0,\\cdots,2m$, og anvend Simpsons formel på hvert delintervall $[x_{2j}, x_{2j+2)}]$. Det resulterer i \n",
"\n",
"\n",
"\\begin{align}\n",
"\\int_a^b f(x)dx &= \\sum_{j=0}^{m-1} \\int_{x_{2j}}^{x_{2j+2}} f(x) dx \n",
"\\approx \\sum_{j=0}^m S(x_{2j},x_{2j+2}) \\\\\n",
"& = \\sum_{j=0}^{m-1} \\frac{h}{3} \n",
"\\left[ f(x_{2j}) + 4 f(x_{2j+1})+ f(x_{2j+2}) \\right] \\\\ &= \n",
"\\frac{h}{3} \\left[ f(x_0) + 4\\sum_{j=0}^{m-1}f(x_{2j+1}) + 2 \\sum_{j=1}^{m-1}f(x_{2j}) + f(x_{2m}) \\right] = S_{m}(a,b). \n",
"\\end{align}\n",
"\n",
"Den sammensatte formelen, her kalt $S_{m}(a,b)$, er kjent som [Simpsons regel](https://wiki.math.ntnu.no/tma4100/tema/numerics?numerisk_integrasjon)\n",
"fra Matematikk 1. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Implementasjon av sammensatt Simpsons metode"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"from numpy import *\n",
"from matplotlib.pyplot import *\n",
"from math import factorial\n",
"newparams = {'figure.figsize': (8.0, 4.0), 'axes.grid': True,\n",
" 'lines.markersize': 8, 'lines.linewidth': 2,\n",
" 'font.size': 14}\n",
"rcParams.update(newparams)"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def simpson(f, a, b, m=10):\n",
"# Finner en tilnærmelse til integralet av funksjonen f\n",
"# over intervallet [a,b] ved bruk av Simpsons metode med m delintervaller.\n",
" n = 2*m\n",
" x_noder = linspace(a, b, n+1) # Jevnt fordelte noder fra a til b\n",
" h = (b-a)/n\n",
" S1 = f(x_noder[0]) + f(x_noder[n])\n",
" S2 = sum(f(x_noder[1:n:2])) # S2 = f(x_1)+f(x_3)+...+f(x_m)\n",
" S3 = sum(f(x_noder[2:n-1:2])) # S3 = f(x_2)+f(x_4)+...+f(x_{m-1})\n",
" S = h*(S1 + 4*S2 + 2*S3)/3\n",
" return S "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Testing\n",
"_Test om koden er korrekt_: Simpsons metode vil integrere 2.gradspolynomer eksakt. f.eks. $f(x)= 4x^3 + x^2+2x-1$. Vis at $\\int_{-1}^2 f(x)dx = 18$. Sjekk om koden gir rett svar."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"I = 18.00000000, S = 18.00000000, feil = 0.000e+00\n"
]
}
],
"source": [
"# Definerer funksjonen\n",
"def f(x):\n",
" return 4*x**3+x**2+2*x-1\n",
"\n",
"I_eksakt = 18.0\n",
"\n",
"# Bruker Simpsons sammensatte metode for å tilnærme integralet\n",
"S = simpson(f, -1, 2, m=1)\n",
"feil = I_eksakt-S\n",
"print('I = {:.8f}, S = {:.8f}, feil = {:.3e}'.format(I_eksakt, S, feil))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"##### Numeriske eksperimenter\n",
"1. Finn en tilnærmelse til $\\int_{0}^{\\pi/2} \\cos(x) dx$. Bruk $m=1,2,4,8$. Finn eksperimentelt hvor stor $m$ må være for at feilen skal bli mindre enn $10^{-6}$. \n",
"2. Gjenta eksperimentet med integralet\n",
"$$ \\int_{0}^8 \\frac{1}{1+16x^2} dx. $$\n",
"(Eksakt løsning er $\\arctan(32)/4$.)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Feilanalyse\n",
"Bruk [Taylorutviklinger](https://wiki.math.ntnu.no/tma4100/tema/taylorpolynomials) rundt midtpunktet $c = (b+a)/2$. La $h=(b-a)/2$. Dette gir\n",
"\n",
"\\begin{align}\n",
"\\int_a^b f(x)dx &= \\int_{-h}^{h} f(c+s)ds =\n",
"\\int_{-h}^h f(c) + sf'(c) + \\frac{1}{2}s^2 f''(c) + \\frac{1}{6} s^3 f'''(c) + \\frac{1}{24}s^4 f^{(4)}(c) + \\dotsm \\\\\n",
"&= 2h f(c) + \\frac{h^3}{3} f''(c) + \\frac{h^5}{60} f^{(4)}(c) + \\dotsm. \n",
"\\end{align}\n",
"\n",
"Hvert andre ledd forsvinner pga. symmetri. \n",
"\n",
"Tilsvarende for kvadraturet $S(a,b)$:\n",
"\n",
"\\begin{align}\n",
"S(a,b) &= \\frac{h}{3}\\left( f(c-h)+4f(c)+f(c+h) \\right) \\\\\n",
" &= \\frac{h}{3}\\left( f(c) - hf'(c) + \\frac{1}{2}h^2 f''(c) - \\frac{1}{6} h^3 f'''(c) + \\frac{1}{24}h^4 f^{(4)}(c) + \\dotsm \\right. \\\\\n",
" &\\qquad + 4f(c) \\\\\n",
" &\\qquad + \\left. f(c) + hf'(c) + \\frac{1}{2}h^2 f''(c) + \\frac{1}{6} h^3 f'''(c) + \\frac{1}{24}h^5 f^{(4)}(c) + \\dotsm \\right) \\\\\n",
" &= 2h f(c) + \\frac{h^3}{3} f''(c) + \\frac{h^5}{36} f^{(4)}(c) + \\dotsm\n",
"\\end{align}\n",
"\n",
"Vi har da følgende uttrykk for feilen: \n",
"\n",
"$$ \n",
" E(a,b) = \\int_a^b f(x) dx - S(a,b) = -\\frac{h^5}{90} f^{(4)}(c) + \\cdots \n",
" = - \\frac{(b-a)^5}{2^5 \\cdot 90} f^{(4)}(c) + \\dotsm\n",
"$$\n",
"\n",
"Dette gir en rimelig ide om størrelsen på feilen, under forusetningen av at $b-a$ er relativt liten slik at høyere ordens ledd kan ignoreres. \n",
"\n",
"Det er mulig, men ikke helt trivielt å vise følgende presise resultat: \n",
"\n",
"##### Setning (feil for Simpsons metode)\n",
"\n",
"_La $f(x) \\in C^{4}[a,b]$. Da fins en $\\xi \\in (a,b)$ slik at_\n",
"\n",
"$$\n",
" E(a,b) = \\int_a^b f(x)dx - \\frac{b-a}{6} \\left[\\;f(a)+4f\\left(\\frac{b+a}{2}\\right)+f(b)\\;\\right] = -\\frac{(b-a)^5}{2880}f^{(4)}(\\xi). \n",
"$$\n",
"\n",
"Dette kan vi bruke for å finne et uttrykk for feilen i sammensatt Simpsons formel. \n",
"\n",
"\\begin{align}\n",
"\\int_a^b f(x)dx - S_{m}(a,b) & = \n",
"\\sum_{j=0}^{m-1} \\left( \\int_{x_{2j}}^{x_{2j+2}} f(x)dx - \\frac{h}{3} \n",
"\\left[ f(x_{2j}) + 4 f(x_{2j+1})+ f(x_{2j+2}) \\right] \\right) \\\\\n",
"& = \\sum_{j=0}^{m-1} -\\frac{(2h)^5}{2880} f^{(4)}(\\xi_j)\n",
"\\end{align}\n",
"\n",
"der $\\xi_{j} \\in (x_{2j}, x_{2j+2})$. Vi kan nå bruke [middelverdisetningen](#gen_middel) til å vise at det fins en $\\xi \\in (a,b)$ slik at \n",
"\n",
"$$\n",
" \\sum_{j=0}^{m-1} f^{(4)}(\\xi_j) = m f^{(4)}(\\xi). \n",
"$$\n",
"\n",
"Ved å bruke $2mh = b-a$ får vi følgende uttrykk for feilen:\n",
"\n",
"$$\n",
" \\int_a^b f(x)dx - S_{m}(a,b) = -\\frac{(b-a)h^4}{180} f^{(4)}(\\xi). \n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Feilestimat\n",
"\n",
"Er det mulig å finne en tilnærmelse til feilen $I(a,b)-S(a,b)$, uten å kjenne $f^{(4)}(x)$? \n",
"\n",
"Anta at intervallet $(a,b)$ er så lite at $f^{(4)}(x)$ er tilnærmet konstant over intervallet. La $H=b-a$. \n",
"\n",
"Da har vi følgende: \n",
"\n",
"\\begin{align}\n",
"I(a,b) - S_1(a,b) &\\approx C H^5, \\\\\n",
"I(a,b) - S_2(a,b) &\\approx 2 C \\left(\\frac{H}{2}\\right)^5. \n",
"\\end{align}\n",
" \n",
"der $I(a,b) = \\int_a^b f(x)dx$, og $S_m(a,b)$ er sammensatt Simpsons formel anvendt på $m$ like store delintervaller av $[a,b]$, og $C = f^{(4)}(\\xi)/2880$, der $\\xi$ er et eller annet punkt i $(a,b)$. Hvilket spiller ingen rolle siden forutsetningen her er at $f^{(4)}(x)$ er nesten konstant.\n",
"\n",
"Ta differensen mellom disse to:\n",
"\n",
"$$ \n",
" S_2(a,b) - S_1(a,b) \\approx \\frac{15}{16}C H^5 \n",
" \\qquad \\Rightarrow \\qquad\n",
" CH^5 \\approx \\frac{16}{15}(S_2(a,b) - S_1(a,b)).\n",
"$$\n",
"\n",
"Sett dette inn i uttrykket over:\n",
"\n",
"\\begin{align}\n",
"E_1(a,b) = I(a,b) - S_1(a,b) &\\approx \\frac{16}{15} (\\,S_2(a,b) - S_1(a,b)\\, ) = \\mathcal{E}_1(a,b), \\\\\n",
"E_2(a,b) = I(a,b) - S_2(a,b) &\\approx \\frac{1}{15} (\\,S_2(a,b) - S_1(a,b)\\, ) = \\mathcal{E}_2(a,b).\n",
"\\end{align}\n",
"\n",
"Siden feilen i $S_2(a,b)$ er omtrent 1/16 av feilen i $S_1(a,b)$, og begge uansett må regnes ut for å finne feilestimatene, så vil vi bruke $S_2(a,b)$ som tilnærmelse. Denne kan gjøres enda bedre ved å legge til feilestimatet. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"###### Eksempel\n",
"Finn en tilnærmelse til integralet $\\int_0^{1} \\cos(x)dx = \\sin(1)$ ved hjelp av Simpsons formel over et og to delintervaller. Beregn feilestimatene $\\mathcal{E}_m$, $m=1,2$ og sammenlign den med virkelig feilen. \n",
"\n",
"Vi får\n",
"\n",
"\\begin{align}\n",
"S_1(0,1) &= \\frac{1}{6} \\big[ \\cos(0.0) + 4\\cos(0.5) + \\cos(1.0) \\big] = 0.8417720923 \\\\\n",
"S_2(0,1) &= \\frac{1}{12} \\big[ \\cos(0.0) + 4 \\cos(0.25) +2 \\cos(0.5) + 4 \\cos(0.75) + \\cos(1.0) \\big] = 0.8414893826\n",
"\\end{align}\n",
"\n",
"Den eksakte feilen og feilestimatet blir\n",
"\n",
"\\begin{align}\n",
"E_1(0,1) &= \\sin(1) - S_1(0,1) = -3.011 \\cdot 10^{-4}, &\n",
"\\mathcal{E}_1(0,1) &= \\frac{16}{15}(S2-S1) = -3.016\\cdot 10^{-4}, \\\\\n",
"E_2(0,1) &= \\sin(1)-S_2(0,1) = -1.840 \\cdot 10^{-5}, &\n",
"\\mathcal{E}_2(0,1) &= \\frac{1}{16} (S2-S1) = -1.885 \\cdot 10^{-5}.\n",
"\\end{align}\n",
"\n",
"I dette tilfellet er det et meget godt samsvar mellom feilestimatet og målt feil. Vi får en enda bedre tilnærmelse ved å bruke\n",
"\n",
"$$ Q = S_{2}(0,1) + \\mathcal{E}_2(0,1) = 0.8414705353607151 $$\n",
"\n",
"som har en feil $\\sin(1)-Q = 4.4945 \\cdot 10^{-7}$. Det gir et mye mer nøyaktig resultat uten ekstra arbeid. \n",
"\n",
"### Implementasjon av Simpsons metode med feilestimat\n",
"Funksjonen `simpson_basis` regner ut en tilnærmelse til integralet \n",
"$$\\int_{a}^b f(x)dx $$\n",
"ved hjelp av Simpsons metode, inkludert et feilestimat. "
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def simpson_basis(f, a, b):\n",
" # Regner ut en tilnærmelse til integralet av f over [a,b].\n",
" # Inn: Funksjonen f, og integrasjonsintervallet [a,b]\n",
" # Ut: S_1(a,b), S_2(a,b) og feilestimatet for S_2(a,b).\n",
" \n",
" # Finner nodene \n",
" c = 0.5*(a+b)\n",
" d = 0.5*(a+c)\n",
" e = 0.5*(c+b)\n",
" \n",
" # Regner ut S1=S_1(a,b), S2=S_2(a,b) og feilestimatet for S2\n",
" H = b-a\n",
" S1 = H*(f(a)+4*f(c)+f(b))/6\n",
" S2 = 0.5*H*(f(a)+4*f(d)+2*f(c)+4*f(e)+f(b))/6\n",
" feilestimat = (S2-S1)/15\n",
" return S1, S2, feilestimat"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Testing\n",
"Test implementasjonen på eksempelet over: \n",
"$$ \\int_0^1 \\cos(x) dx. $$ \n",
"og kontroller at svarene blir de samme som tidligere. "
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"S1 = 0.84177209, S2 = 0.84148938\n",
"Feil i S2 = -1.840e-05, Feilestimat for S2 = -1.885e-05\n"
]
}
],
"source": [
"def f(x):\n",
" return cos(x)\n",
"\n",
"a, b = 0, 1\n",
"I_eksakt = sin(1) # Eksakt løsning\n",
"\n",
"S1, S2, feilestimat = simpson_basis(f, a, b)\n",
"\n",
"print('S1 = {:.8f}, S2 = {:.8f}'.format(S1, S2))\n",
"print('Feil i S2 = {:.3e}, Feilestimat for S2 = {:.3e}'.format(I_eksakt-S2, feilestimat))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"##### Numerisk eksperiment\n",
"Gitt det ubestemte integralet (med ubestemt løsning):\n",
" $$ \\int \\frac{1}{1+16x^2} dx = \\frac{\\arctan(4x)}{4}. $$ \n",
"1. Bruk `simpson_basis` for å finne en tilnærmelse til integralet over intervallet $[0,8]$. Skriv ut $S_2(0,8)$, feilestimatet $\\mathcal{E}_2(0,8)$ og målt feil $E_2(0,8)$. Er feilestimatet pålitelig? \n",
"2. Gjenta eksperimentet over intervallene $[0,1]$ og $[4, 8]$. Legg merke til forskjellen i målt feil over de to intervallene.\n",
"3. Gjenta eksperimentet igjen over intervallet $[0,0.1]$"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"S1 = 0.09513705, S2 = 0.09512722\n",
"Feil i S2 = -6.282e-07, Feilestimat for S2 = -6.550e-07\n"
]
}
],
"source": [
"def runge(x):\n",
" return 1/(1+16*x**2)\n",
"\n",
"a, b = 0, 0.1\n",
"I_eksakt = 0.25*(arctan(4*b)-arctan(4*a)) # Eksakt løsning\n",
"\n",
"S1, S2, feilestimat = simpson_basis(runge, a, b)\n",
"\n",
"print('S1 = {:.8f}, S2 = {:.8f}'.format(S1, S2))\n",
"print('Feil i S2 = {:.3e}, Feilestimat for S2 = {:.3e}'.format(\n",
" I_eksakt-S2, feilestimat))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Observasjoner av eksperimentet:\n",
"1. Over intervallet $[0,8]$: Stor feil, dårlig feilestimat.\n",
"2. Over intervallet $[0,1]$: Stor feil, dårlig feilestimat. Over intervallet $[4,8]$: Liten feil, greit feilestimat. \n",
"3. Over intervallet $[0,0.1]$, liten feil, greit feilestimat.\n",
"\n",
"_Forklaring_\n",
"\n",
"Fra [setningen om feil i Simpsons metode](#feil_simpson) har vi\n",
"$$\n",
" E(a,b) = -\\frac{(b-a)^5}{2880}f^{(4)}(\\xi). \\tag{1}\n",
"$$\n",
"\n",
"Så vi må forvente at feilen er stor når $f^{(4)}(x)$ er stor. Videre antok vi i utledningen av feilestimatet at $f^{(4)}(x)$ var nogenlunde konstant over intervallet. \n",
"\n",
"Så det kan være en god grunn til å ta en titt på den 4. deriverte av vår funksjon: \n",
"\n",
"$$ f(x)=\\frac{1}{1+16x^2} \\qquad \\Rightarrow \\qquad \n",
" f^{(4)}(x) = 6144 \\frac{1280 x^4 - 160x^2 +1}{(1-16x^2)^5} $$\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAgYAAAEiCAYAAACC6hodAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3XmcZGV97/HPt7unZ4cBhpkBFAZk\nERFEHG6AsDTGcbmaRL3eQEQFk4AG0UtQDGjESaIGNxa3C5gFRQlEzDWCymKgkbA6RJaRVZxhmw0Y\nZqCHnulZfveP51R1caZ6P9XVp/y+X696Vdc5zznn91RXd/3Os5yjiMDMzMwMoK3ZAZiZmdnE4cTA\nzMzMqpwYmJmZWZUTAzMzM6tyYmBmZmZVTgzMzMysyomBFUpSt6TucTzepZKWjcNxuiSFpK4y7Nda\nj6SFku6W1Jt9ZuY36DgnZfs/rIB9jev/AyuGEwMblpp/FvUeFzU7PhsdSa+QtEjSwc2OpVYWU+1n\nbJOkZZK+JmlWs+Mbb5K2B/4NEPBR4P3AM00NylpWR7MDsNJZBDyWW/ZIzc9vHr9QWsIvgKlAX5OO\n/wrgs8Ay4J4mxTCY04B1wHTgD0hfim+QdGT8bl2d7fXALODvIuJHzQ5mBPz/oIScGNhIXRcRdwy0\nMiIK/YKTNC0iXipynxNBpV4RsRXY0ITjdwJbxvu4o/DDiFiZ/XyxpCuBPwEOBe5qXljjbk72vLap\nUYxQ0f8PbHy4K8EKVa9PUclpku6TtEHSc5K+L+kVdbZ9SNJBkm6UtB74Vs36D0p6JNvHvZL+cIAY\nhnW8Qeqwn6TrJL0kaaWkLwOdA5RdIOkaSWuzvt+7JL0jV6bSDXOspAskrQTWZ+teNsZA0iey16+q\nc6xt1knaR9IVkp6RtFHS/ZI+mNuucowTJH1G0hNAL6k5+vas2L/UNNsvGuS96ZT0d5IWS3o+q/Od\nkv4oV26JpFsG2MeA64bhF9lz7XswP4v7pDrHWibp0prXld/FUZK+IGlFVofrJe1ZZ/tTJP0m+xzd\nI+kdqjOuZQSf8b0l/Vt23A2SnpT0A0m7DVTh7O/pyuzlTVn83ZV1+b+3bPkiSZFbFpIukvQ2Sb/K\njv+opD8d6Ng1287IjvWspEOyZXMl/WNWh42Slmd/CwfWxp6PT9JUSedmv5uNkpZK+pykyblyyyRd\nK+l/SPqv7Pf0pKTTh4rXxsYtBjZS20uanVv23BDNuv8X+AvgMtIX/TxSk/ARkl4fEbVnQbOA64F/\nB64gO0OSdCLwz8AvgW8Cc7P9PTnG472MpDnAzcA04KvAc8AHgIV1yh4DXAfcD3wO2AgcB/xY0nsi\n4t9zm3w9q88XgO0HCOFK4EvA8cDnc+uOBxZHxGPZ8fcHbgWeBc7L9v124J8lbR8RF+S2PxsI4EJS\nX/XPSF1Di4BLgMqX9X0DxAawHfDhLM5/JnWDnAD8SNLbIuK6rNwVwN9JekVEPFXZWNIBwAGkLoLR\nmJ89rxnl9hXnk5KjfwBmA58Avg8cUSkg6RTgYuAO0u9uDukz9VR+ZwzjMydpEunzMpX0GV4B7AK8\nhdSl8/QAsX4eWAJ8hPTZeRBYNcp6Hwa8E7gI+Kcs5u9JuiciHqy3gdKYjmtJ731XRCzJVl0FHAh8\nA1gK7AwcDexH+puoty+R/rbfCnwXuBM4Evg06XPxrtwmewJXA5eSfj/HAedLeiAirh9Z1W3YIsIP\nP4Z8ACeRvlTqPWbXlOsGumteH5GVOSm3v4OATcBnctsG8LFc2Q5gJfBrYErN8jdn5ZeN5ngD1POr\n2fZH1SybBjyaLe/Klgl4CLgRaKsp20Y6C3+sznt3J9CRO15X7X6zZbcA9+fKvSor9/GaZdcDDwDT\ncmWvBF4ApueO8URlWU3Zw+q9X4O8P+3A5NyyziyOGwaLN1v+OWAzMHeI4yzKtn8N6Yt7D+CDwEuk\nL9Taz8H8gepAGjtxaZ3fxc2539vp2fIDsteTSIP77gE6a8r9wWg/c8DrsnLvGcXf3/H5z0m9v7f8\n+5dbFlk8+9csm0tKaL9c5z06LHvv/5uUgO9bU2b7rMwnhoj7ZfEB78i2+/tcufOz5W/K/e7yyyaT\nkqIfjPQ99GP4D3cl2Eh9jHT2XPtYN0j5PwF6gJ9Kml15AMtJX7ZvzJXfTDp7rXUo6R/YRRFR7Y+P\ndMbwwBiPl/cO4L8jotrUHWmMQz6m15HOjL4H7FhznB1JZ+J7Sdojt823I2LzEMeHdLb9WkmvqVlW\n+WL4NwBJOwBvyl5Py9X1Z8BMYEFuv9+NiPXDOP6AImJLRGzMYuiUtCOpFeFm4A015R4DFpPO8God\nB9wUEcM94/016Qt6GamF4iHgbbWfg1G6ONL4joqbs+e9sudDSV+Kl0RNP3lE/GcWU63hfuZeyJ7f\nImn6GOMfrZuipmUg+z08RH+9a+1Cel9mAUdHRO0g4w2kJKMr+wwMVyUx+Gpu+Zey57fnlj8aET+v\niXcjqQWnXrxWEHcl2Ej9MgYZfFjHvsAMBm76zHdBLK/zT7/yBftwne0fBg4Zw/Hy9iA1ddY7Tq19\ns+d/GmRfc4DHa17nZ3MM5Aek5v7jgXOyZccDt0ZEpetkH1KrxWezx0DHrzXc4w9K0l8AfwXsn8VQ\nkX9vrwC+ImmviPitpDcAewPnjuBwfwI8T6rLR0lNyxtHG3uNx3Ovn8+eK19ylc/co3W2fYRRfOYi\nYqmkC4H/A7xP0q3ANcD3IuLZkYU/avl6Q6p7vS/372bPr6n53AHpC1rSWaQv9FWS7gR+SqrLE4Mc\nfz6wMnLdeRGxQtJa+ruKhor3oEGOYWPkxMAarY3UT3/8AOvzZ7C9dcpUvnzqfakr93qkx6tnuMcB\nOAu4e4D95JOJenXb9uARqyXdSDq7PidrOXgtL++Xrxz/fNI/5Hry/bzDOv5gsoFq3wZ+DHwRWE1q\n5fkg8N5c8SuBL5N+F1/InjdRP/EayC2RzUqQ9B+k8Q+XS3pDzRn/YMle+wDLB5qRkf89D6fMsD9z\nEXG6pG8Df0TqCvsK8DeSjomIfEvEcESdeKCYel9F6lb4BCmZefmBI86T9P+APya1HH4G+JSkP85a\nVkaqXgxj+T3ZKDkxsEZ7jPRP486IeHGU+1iWPb8auCG3bt/c67Ee7/HsOHn1jgPwYm1TZ4GuAP4p\nGwH+LtI/yKtq1v82e948xuOP9FoAx2XHfmdEVLdVbiYEQEQ8Jem/gOMl/QPp7P+6iHg+X3ZYgUas\nl/RZ0gC/44HLs1WVgYgvu/BRNsp9l9Eci/4z1X3Y9jO3T+71iD5zWQLwa+AfJB1ESixPB04eRZzP\nU79Zff4o9pV3MSnOL0taHxGfyheIiKXABcAFkl4J/IqULA+UGCwD3ixpVm2rgaR5pHELywqI28bI\nYwys0a4gfc4W5VdkU7zyMxzqWUw6M/2QpCk127+ZNDityOP9BDhE0lE1200DTsmVu5vUzPwJpavS\n5Y+18xDHGcq/ky56dBx1+uUjYjVp4OMp2T/k0R6/cja7wzDLV87gqmdskvZi29HkFf9KGrn+58Du\npN/PWFxBGkR5djbCnezL+Bng2FzZDzPwmfNQFpNme5yidM0HACT9AWn0fD6mIT9zkraTlD8Ze5DU\nkjPc9z/vN8D+kubWHHM30syDMYuIrwB/S3q/q4mBpGmSpubKPkn6Ox2sLteQPjt/lVt+Zvb8kzEH\nbWPmFgNrqIi4RdLXgDMkvY40MO4lUl/xu0h99IP2OUfEJklnZ2V/Ien7pD7nj5DOaGYUeLwvAu8D\nrsn2U5mu+LIuiIjYKunPSNPPHpD0z6SzzF2Aw0l91PkvkGGLNL3tWuAvSQMJv1in2F+SpivelzVP\nP0oaMHcIaQrcdsM41KOkQXF/KakHeBFYEv1T0vJ+DLybNCXzx8BuwKmkfvfX1Sl/FfA10nTKXuA/\nhhHTgCJic9ZP/1XSQLars1UXk5rk/4U0OG0BaQbBqPruI6JP0mdI0xBvlvSv9H/mljC6z9wbgW9K\nuorUzSRS0jeT0SdM/wh8HLg++wzMIn0u8uMgRi0iFmWDJT8vqScivkZqQbtR0g9If4Mbgf9JGndy\n5sB74yekqY/nSNqdNP34CNKU1/9oUOubjVSzp0X4UY4HNVOYhijXTf3pUx8g/cNeT/ryeZA0/3m/\n3LYPDbLvPyd9kW0k9TX/IWl+87LRHG+Q4+xPaj7uJQ0o+wr9UyO7cmUPIPWlr87iepL0ZfWemjID\nvnfUma5Ys64yE6EP2GGAWHcnfTk8nZVbDvwcOKXOMY4fYB9/SBqP0JeVWzTE+3MmqTthA+lL8n3U\nmR5XU/7abL/DnmJG/3TFeXXWzSRds+G2mmVTSF/ia7Lf+U9ITezLqD9d8bDcPudTf8rhh0ldBRtI\nUxffQUp2HhzpZ46UKPxj9hl+KYv1FuCPh/F+1J2umK37E1Ii0Jcd80/r/T6y7S8a6m92kPfoW8BW\n4M+AnUjXdnggq+sLpFaWPxvq/wHpOg5fJLX89GW/o8+x7TTYZcC1deK9lDp/834U91D2RpsVQumK\ndhsj4k3NjsWsESTdC6yOiG0uemUv5/8H5eQxBla0XRhl863ZRCJpSmUcQ82yPyBNlbuxOVGVjv8f\nlJDHGFghJB1B6nt+FfX7w83K5jDg69mYgBWkbqMPkbptfKvxQfj/Qbk5MbCinAy8jTR16V+aHItZ\nEZaRxlJ8mNSnvg74f8DZMcopl79D/P+gxDzGwMzMzKp+J1sMZs+eHfPnzy9sf+vXr2f69GZd+rxY\nrsvE1Cp1aZV6gOsyEbVKPaD4utx9993PRsSwrm/yO5kYzJ8/n8WLFxe2v+7ubrq6ugrbXzO5LhNT\nq9SlVeoBrstE1Cr1gOLrIqnefSfq8qwEMzMzq3JiYGZmZlVODMzMzKzKiYGZmZlVOTEwMzOzKicG\nZmZmVuXEwMzMzKqcGIzRN258lM/e1sviZWuaHYqZmdmYNSQxkLSLpO9IekbSBkkPSDqmZr0kLZK0\nXFKvpG5JB+T2sYOkyyStyx6XSZqVK3OgpJuzfTwt6Zz83dAa7ck1vTz+wlbue2rdeB7WzMysIQpP\nDLIv71sBAW8H9gc+CqyuKfZJ4OPZ8kOzdTdImllT5nLgENKNON6a/XxZzXG2A24AVmX7+BhwJnBG\n0XUazAG7bQfAr5e/MJ6HNTMza4hGXBL5k8CKiPhAzbKllR+yM/rTgXMj4ofZshNJycF7gYsl7U9K\nBo6MiNuyMh8CbpG0X0Q8DJwATANOjIheYEm23RmSzotxujvUfnNTLvPYMz3jcTgzM7OGakRXwjuB\nOyVdKWm1pHsknVbTxL8nMA+4vrJB9sX+C+CIbNHhQA9wW81+bwXW58rckm1bcR2wKzC/2CoNbKcZ\nnQC80LtpvA5pZmbWMI1oMdgLOBU4HzgXOBj4erbuG6SkAFIXQK1VwG7Zz/OAZ2rP+iMiJK2u2X4e\n8FSdfVTWLa1dIekU4BSAuXPn0t3dPdJ61bV241YAnn1hfWH7bKaenp6WqAe4LhNRq9QDXJeJqFXq\nAc2tSyMSgzZgcUScnb3+laR9gI+QEoOKfFO/csvqdQUMVUYDLCciLgEuAViwYEEUddeqDZu2wE3X\n0rtFHHPMMYzz2MfC+e5kE1Or1KVV6gGuy0TUKvWA5talEV0JK4AHcsseBHbPfl6ZPc/LlZlD/xn/\nSmBO7QyD7Oedc2Xq7QO2bY1omCmT2pnUBpu2BBs2bR2vw5qZmTVEIxKDW4H9csv2BSr3gl5K+lJf\nWFkpaQpwFP1jCm4HZpDGEVQcDkzPlTkq27ZiIbAcWDbWSozE9Ekpf1nncQZmZlZyjUgMzgcOk/Rp\nSXtL+t+kqYTfhDRWALgAOEvSuyW9FriUNNjw8qzMg8C1pBkKh0k6HLgYuCabkUBW9iXgUkmvlfRu\n4Cxg3GYkVEzLOmScGJiZWdkVPsYgIn4p6Z3AF4DPAE9kz9+qKfYlYCopWdgBuBN4c0S8WFPmBOBr\n9M9e+DFwWs1x1klamO1jMfA88FXgvKLrNJTJHWnoQ++mLeN9aDMzs0I1YvAhEfET4CeDrA9gUfYY\nqMwa4H1DHOd+4OhRBVmgzqzd5aW+zc0NxMzMbIx8r4QCTG5PYww2uMXAzMxKzolBATrb03Nvn2cl\nmJlZuTkxKEClxcBdCWZmVnZODApQaTFwV4KZmZWdE4MCTK50JTgxMDOzknNiUIDOaleCEwMzMys3\nJwYF6HSLgZmZtQgnBgWY3JZaDHrdYmBmZiXnxKAAk7IWg77Nnq5oZmbl5sSgAB3Zu+jEwMzMys6J\nQQE6sq6EjVucGJiZWbk5MSjApOxd3LjJiYGZmZWbE4MCVBKDPrcYmJlZyTkxKEClK6Fvs2clmJlZ\nuTkxKMAkDz40M7MW4cSgAB3uSjAzsxbhxKAA/V0JTgzMzKzcnBgUoDorwYmBmZmVnBODAvgCR2Zm\n1iqcGBRgkrsSzMysRTgxKIBbDMzMrFU4MShAdYyBZyWYmVnJOTEoQG2LQUQ0NxgzM7MxcGJQgDap\nf8qiWw3MzKzEnBgUpDNrNvA4AzMzKzMnBgWZ7MTAzMxagBODglRbDNyVYGZmJebEoCDuSjAzs1bg\nxKAgne1ODMzMrPycGBSks6Md8P0SzMys3JwYFKTSleDEwMzMysyJQUE8K8HMzFqBE4OCTPasBDMz\nawFODAriwYdmZtYKGp4YSPqUpJD0jZplkrRI0nJJvZK6JR2Q224HSZdJWpc9LpM0K1fmQEk3Z/t4\nWtI5ktToOtXj6YpmZtYKGpoYSDoMOBm4L7fqk8DHgY8ChwKrgRskzawpczlwCPA24K3Zz5fV7Hs7\n4AZgVbaPjwFnAmc0oi5D6b/A0ZZmHN7MzKwQDUsMJG0PfB/4c+D5muUCTgfOjYgfRsQS4ERgJvDe\nrMz+pGTglIi4LSJuBz4EvEPSftmuTgCmASdGxJKI+CHwReCMZrQaVLoSNm5yi4GZmZVXI1sMLgGu\niogbc8v3BOYB11cWREQv8AvgiGzR4UAPcFvNdrcC63Nlbsm2rbgO2BWYX0wVhq/SYrBpq2+7bGZm\n5dXRiJ1KOhnYG3h/ndXzsudVueWrgN1qyjwTEdVv2YgISatrtp8HPFVnH5V1S3MxnQKcAjB37ly6\nu7uHW50h9fT0sHrlRgAefOgRujcsHWKLiaunp6fQ96aZXJeJp1XqAa7LRNQq9YDm1qXwxCBr6v8C\ncFRE9A1SNH9qrdyyeqfeQ5XRAMuJiEtIrRgsWLAgurq6BgltZLq7u9lzj7nw+G/ZY8+96DrmVYXt\ne7x1d3dT5HvTTK7LxNMq9QDXZSJqlXpAc+vSiK6Ew4HZwBJJmyVtBo4BTs1+fi4rNy+33Rz6z/hX\nAnNqxwpkP++cK1NvH7Bta0TDTWpPoW7ydQzMzKzEGpEY/Ag4EDi45rEYuCL7+RHSl/rCygaSpgBH\n0T+m4HZgBinJqDgcmJ4rc1S2bcVCYDmwrMgKDcekynUMtniMgZmZlVfhXQkRsRZYW7tM0npgTTYD\nAUkXAJ+W9BApUfgb0mDDy7N9PCjpWuDibLyCgIuBayLi4Wy3lwOfBS6V9DlgX+As4G9rxyaMl0pi\n4BYDMzMrs4YMPhyGLwFTgW8COwB3Am+OiBdrypwAfI3+2Qs/Bk6rrIyIdZIWZvtYTJoS+VXgvIZH\nX0dluuImX+DIzMxKbFwSg4joyr0OYFH2GGibNcD7htjv/cDRYw6wAB5jYGZmrcD3SijIpA6PMTAz\ns/JzYlAQjzEwM7NW4MSgIJ1ODMzMrAU4MSiIWwzMzKwVODEoSGXwYd9mjzEwM7PycmJQkMrgQ7cY\nmJlZmTkxKIjHGJiZWStwYlAQjzEwM7NW4MSgINUxBr6OgZmZlZgTg4JM8iWRzcysBTgxKEinBx+a\nmVkLcGJQEI8xMDOzVuDEoCD9N1HyGAMzMysvJwYFqUxX7HOLgZmZlZgTg4K4K8HMzFqBE4OCVK98\n6FkJZmZWYk4MCuIxBmZm1gqcGBRkUlv/GIMIJwdmZlZOTgwK0tYmOtpSq8HmrU4MzMysnJwYFMgD\nEM3MrOycGBSoOs5gs1sMzMysnJwYFKhyWWRfy8DMzMrKiUGB3JVgZmZl58SgQE4MzMys7JwYFKj/\nWgZODMzMrJycGBSo0mLQ58GHZmZWUk4MClQZfOgWAzMzKysnBgXyGAMzMys7JwYFqowx8HRFMzMr\nKycGBepvMfAYAzMzKycnBgXqbPetl83MrNycGBTIYwzMzKzsnBgUaJIviWxmZiXnxKBA/Rc48hgD\nMzMrp8ITA0lnS/qlpBckPSPpakmvzZWRpEWSlkvqldQt6YBcmR0kXSZpXfa4TNKsXJkDJd2c7eNp\nSedIUtF1Gq5OdyWYmVnJNaLFoAv4FnAE8EZgM/BzSTvWlPkk8HHgo8ChwGrgBkkza8pcDhwCvA14\na/bzZZWVkrYDbgBWZfv4GHAmcEYD6jQsHmNgZmZl11H0DiPiLbWvJb0fWAf8PnB1dkZ/OnBuRPww\nK3MiKTl4L3CxpP1JycCREXFbVuZDwC2S9ouIh4ETgGnAiRHRCyzJtjtD0nkRMe7t+f2XRHZiYGZm\n5TQeYwxmZsd5Pnu9JzAPuL5SIPti/wWplQHgcKAHuK1mP7cC63Nlbsm2rbgO2BWYX2gNhmlSh8cY\nmJlZuRXeYlDHhcA9wO3Z63nZ86pcuVXAbjVlnqk964+IkLS6Zvt5wFN19lFZt7R2haRTgFMA5s6d\nS3d392jqUldPTw/d3d0sf6oPgEd+8xjdPFnY/sdTpS6twHWZeFqlHuC6TEStUg9obl0amhhIOg84\nktQlsCW3On9ardyyeqfdQ5XRAMuJiEuASwAWLFgQXV1dg8Y+Et3d3XR1dXHflkfhsUd4xe570NW1\nX2H7H0+VurQC12XiaZV6gOsyEbVKPaC5dWlYV4Kk84E/Bd4YEb+tWbUye56X22QO/Wf8K4E5tTMM\nsp93zpWptw/YtjViXFTHGHjwoZmZlVRDEgNJF5IGEr4xIh7KrV5K+lJfWFN+CnAU/WMKbgdmkMYR\nVBwOTM+VOSrbtmIhsBxYVkhFRqh6HYPNHmNgZmbl1IjrGHwT+CCpteB5SfOyxwxIYwWAC4CzJL07\nu8bBpaTBhpdnZR4EriXNUDhM0uHAxcA12YwEsrIvAZdKeq2kdwNnAU2ZkQDQ2eHpimZmVm6NGGNw\navb8n7nlfwssyn7+EjAV+CawA3An8OaIeLGm/AnA1+ifvfBj4LTKyohYJ2lhto/FpFkPXwXOK6oi\nI1XpSti81YmBmZmVUyOuYzDklQezM/pF9CcK9cqsAd43xH7uB44eWYSN038dA3clmJlZOfleCQXq\nv1eCWwzMzKycnBgUyPdKMDOzsnNiUCDfK8HMzMrOiUGBJnVUrmPgMQZmZlZOTgwK1H8dA7cYmJlZ\nOTkxKJDHGJiZWdk5MSiQxxiYmVnZOTEoUP+9EjzGwMzMysmJQYE6O3wdAzMzKzcnBgVyV4KZmZWd\nE4MCVRMDz0owM7OScmJQII8xMDOzsnNiUCBPVzQzs7JzYlCgSR58aGZmJefEoEAefGhmZmXnxKBA\nHW2VFoMgwuMMzMysfJwYFEhSzTgDJwZmZlY+TgwKVr2RkrsTzMyshJwYFKxy62UnBmZmVkZODArW\n0Va5loETAzMzKx8nBgXrbO8fgGhmZlY2TgwKVu1K8GWRzcyshJwYFMzXMjAzszJzYlCw/vslODEw\nM7PycWJQMI8xMDOzMnNiUDB3JZiZWZk5MShYNTHw4EMzMyshJwYFq8xK8BgDMzMrIycGBfMYAzMz\nKzMnBgXzGAMzMyszJwYFc2JgZmZl5sSgYNXrGHjwoZmZlZATg4J1dniMgZmZlZcTg4K5K8HMzMrM\niUHBnBiYmVmZdTQ7gLGSdCpwJrAL8Gvg9Ii4pVnxjOVeCVu2Btfct5yf3r+CR1f18FLfFqZPbmfe\n9lPYc/Z09pw9g71mT2f3naax26ypTJnUXnT4Zmb2O67UiYGk44ALgVOB/8qefybpNRHxRDNiql7H\nYPPIxhisWd/Hyd9dzN2PP7/NuseeWc+tv3nuZcskmLfdFF65wzReseNUZs+YzI7TO9lxWic7Tu9k\n1rRJTO1sZ1pnB1MntadHZzuT2oWk0VfQzMxaWqkTA+AM4NKI+Hb2+qOS3gr8JXB2MwIaTVfCxs1b\nOOlf7uK+p9Yxd7vJnNq1N4fttRMzpnTQs2Ezy9f28ttn17P02R6WPrueJ9f08vTaXlas28CKdRu4\na9nw42tvE5PaRUdbG+1toqNNtGXP7W1i08YNzLi7m462Nir5gyRESkaqy1B6XdlxbRnqbEd1w+qy\nRlu7tpeLHrm98P2K8U+s1q7t5eJH7hjzfpqdEz7/fC+XPDq6ejQ79rznn+/l278Z++9kImiVurRK\nPSDVZc8D17PHTtPH/dilTQwkdQJvAL6SW3U9cESd8qcApwDMnTuX7u7uwmLp6emp7u/Jx/sAeGzZ\n43R3rxjW9lc90sd9T21i9lRx1iFt7NC3jBUPLeuPHXgV8KpZwKy0ZMvWqazZEDzTGzzbu5UX+yJ7\nwIubgvV9Qd9W6NsSbNzS/7xla7BlawADJy6rXlo/indhglqzptkRFGfNc0OXKYNWqQfAc67LhNMq\n9QBuvvUOdt9u/LuMS5sYALOBdmBVbvkq4E35whFxCXAJwIIFC6Krq6uwQLq7u6ns77GOpfDIA8zb\ndTe6ug4YctvVL2zg5/95EwAXnXgYC+bvWFhc9WzaspUtW4PNW4MtW4LNW2tebw1uvf0OFhx6KJu2\nBBEQpGeg/5nKusryIKrro1q2siyyDauvGZ+pnPfecy+vO/h1xe60SbNQ77n3Xg5+3djqMhEm0N57\n7728bhT1iIkQfM5o6zIRtUpdWqUekOry7rccw4zJ4/81XebEoCL/L0N1lo2b/nslDK8r4bI7HmfD\npq28+TVzG54UQOrqGGzM4rzpbew9Z2bD4xgPfU+2c8SrZjc7jEL0PdXOEXuXvy6bnmrn91ugHgCb\nn27nyH1cl4mkVeoBqS7NSAqg3NMVnwW2APNyy+ewbSvCuOm/7fLQucmWrcFVdz8FwJ8duWdD4zIz\nMxuO0iYGEdEH3A0szK1aCNxeg8JmAAAPK0lEQVQ2/hElIxl8eNfSNaxYt4E9dprG7+3Z+NYCMzOz\noZS9K+E84DJJdwG3Ah8GdgUualZAkzqGfx2Dmx5eDcBbD5jnKYRmZjYhlDoxiIgrJe0E/A3pAkdL\ngP8ZEY83K6aRjDG46aGUGByz384NjcnMzGy4Sp0YAETEt4BvNTuOiv6uhMHHGKx+cQOPru5hWmc7\nC/ZwN4KZmU0MpR1jMFENd4zBr55YC8DBr5xFZ4d/DWZmNjH4G6lg1XslbB5eYvD63Wc1PCYzM7Ph\ncmJQsM6O4Y0x+O8n0j0RDtl9h4bHZGZmNlxODAo2nDEGEcGDy18A4KBXuMXAzMwmDicGBRtOV8Ly\ndRt4ceNmdpreyc4zJ49XaGZmZkNyYlCwycO4jsHDK1NrwX7zWuPSw2Zm1jqcGBRscnYjgg2btgxY\n5qGVLwKw71wnBmZmNrE4MSjYlKzFYOMgXQmPruoB4NVuMTAzswnGiUHBpgyjxWDps+sB2HP29HGJ\nyczMbLicGBSsMsZgw6YtxAA3kX9yzUsA7LGTEwMzM5tYnBgUrKO9jY42sTXqT1ns2biZ59b30dnR\nxhzPSDAzswnGiUEDVLoTNm7etjvh8edSN8LuO06jrc13VDQzs4nFiUED9HcnbDsA8Ynnsm6EHaeN\na0xmZmbD4cSgAQYbgPh4Nr5g952cGJiZ2cTjxKABJk+qTFms15XgFgMzM5u4nBg0wJSOSotBna6E\nNWmMgWckmJnZROTEoAEGazF4IutKeKVbDMzMbAJyYtAAA7UYbN0arFq3EYDdZk0d97jMzMyG4sSg\nAaZM6r/IUa3n1vfRt2Urs6ZNYmpnezNCMzMzG5QTgwaY3FG5jsHLWwxWrtsAwLztpox7TGZmZsPh\nxKABBmoxWLGuF4BdtndiYGZmE5MTgwbov45BrsXghdRisIvHF5iZ2QTlxKABBrrA0fK1WWLgrgQz\nM5ugnBg0QOWSyNuOMUhdCfPclWBmZhOUE4MGmDxAi8GKbPDhru5KMDOzCcqJQQNUBx/mLnBUGWPg\nFgMzM5uonBg0QHW6Ys3gw4iothh4VoKZmU1UTgwaYEqdSyKvWd9H3+atbDelg2mdHc0KzczMbFBO\nDBqgcknk3r7+xKA6VXF7jy8wM7OJy4lBA0yfnFoE1tckBquyxGCuuxHMzGwCc2LQADOnpMSgZ8Pm\n6rKV2c2T5m03uSkxmZmZDYcTgwaYkbUY9GysTQyyaxj44kZmZjaBOTFogBlT6iQG7kowM7MScGLQ\nADOzFoMXa7sSXqh0JTgxMDOziavQxEDSjpK+LukhSb2SnpT0fyXtlCu3g6TLJK3LHpdJmpUrc6Ck\nm7P9PC3pHEnKlflfkh6QtDF7fleR9Rmt/haDTdVlq7JrGMx1YmBmZhNY0S0GuwK7AZ8EDgTeBxwN\n/Guu3OXAIcDbgLdmP19WWSlpO+AGYBVwKPAx4EzgjJoyhwNXAt8HDs6efyDp9wqu04hNndROe5vY\nsGkrm7akixz5qodmZlYGhV5pJyKWAO+uWfQbSWcC10jaLiJekLQ/KRk4MiJuA5D0IeAWSftFxMPA\nCcA04MSI6AWWZNudIem8iAjgdOCmiPh8dqzPSzo2W/6nRdZrpCQxY3IH63o3sX7jZqZMamdd7yYm\ntYsdp3U2MzQzM7NBjccYg+2AjcBL2evDgR7gtpoytwLrgSNqytySJQUV15FaJObXlLk+d6zravbR\nVDNqxhmszLoR5sycQlubBtvMzMysqRp6bd5s3MDfA9+OiMpIvHnAM9lZPwAREZJWZ+sqZZ7K7W5V\nzbql2fOqOmXmUYekU4BTAObOnUt3d/doqlRXT0/PNvvT5pQMdN96B+s3papOZWOhx22EenUpK9dl\n4mmVeoDrMhG1Sj2guXUZVmIg6XPAp4codmxEdNdsMx24GniaNOagVrAt5Zbny6jO8npl6u2biLgE\nuARgwYIF0dXVVa/YqHR3d5Pf37wHb+Opnud59YEHs3xtL9x1D/u9ci5dXYcUdtxGqFeXsnJdJp5W\nqQe4LhNRq9QDmluX4bYYXAB8b4gyT1R+kDQD+Gn28h0RsaGm3EpgjiRVWg2y2QY7098CsJJtz/zn\nZM9Dlcm3IjTFjJqrH670jAQzMyuJYSUGEfEs8OxwykqaCfyMdPb+1ojoyRW5HZhBGiNQGWdwODC9\n5vXtwBclTalJKhYCy4FlNWUWAl+u2fdCXj52oWmqYww2bq65gZITAzMzm9iKvo7BTNKAwB2Ak4Dp\nkuZlj06AiHgQuBa4WNJh2bTDi4FrshkJkKYzvgRcKum1kt4NnAWcVzM24ULgjZLOlvRqSWcDx5Ja\nN5qu9n4JvoGSmZmVRdGzEt4AHAa8BngEWFHzqJ0tcAJwLymJuC77+f2VlRGxjnT2vyuwGPgm8FXg\nvJoytwHHAycC9wEfAI6LiDsLrtOozJwyCYB1vZt4ck2aXLHbLCcGZmY2sRV9HYNu+gcJDlZuDeni\nR4OVuZ90caTBylwFXDWCEMfN7BnpegXP9mxk2bPrAZi/0/RmhmRmZjYk3yuhQXaemW6v/PDKF3lx\n42ZmTulgx+m+uJGZmU1sTgwaZOcZqdvgrmVrgNRakLvVg5mZ2YTjxKBB5myXWgz6Nqd7Jeyx07Rm\nhmNmZjYsTgwa5JU7TKO2gWDP2R5fYGZmE58TgwaZ2tnO7jv2txIc9IpZg5Q2MzObGJwYNFBtMnDo\n/B2aGImZmdnwNPQmSr/rzli4L4+uepETDtuDWb7dspmZlYATgwbac/Z0rj190EsxmJmZTSjuSjAz\nM7MqJwZmZmZW5cTAzMzMqpwYmJmZWZUTAzMzM6tyYmBmZmZVTgzMzMysyomBmZmZVSkimh3DuJP0\nDPB4gbucDTxb4P6ayXWZmFqlLq1SD3BdJqJWqQcUX5c9ImLn4RT8nUwMiiZpcUQsaHYcRXBdJqZW\nqUur1ANcl4moVeoBza2LuxLMzMysyomBmZmZVTkxKMYlzQ6gQK7LxNQqdWmVeoDrMhG1Sj2giXXx\nGAMzMzOrcouBmZmZVTkxMDMzsyonBmZmZlblxGCMJJ0qaamkDZLulnRUs2MaKUlHS/qxpKclhaST\nmh3TaEg6W9IvJb0g6RlJV0t6bbPjGg1JH5F0X1aXFyTdLuntzY5rrCR9KvuMfaPZsYyGpEVZ/LWP\nlc2OazQk7SLpO9nfygZJD0g6ptlxjZSkZXV+JyHpJ82ObaQktUv6+5rvlKWSPiepYzzjcGIwBpKO\nAy4EvgC8HrgN+Jmk3Zsa2MjNAJYA/wfobXIsY9EFfAs4AngjsBn4uaQdmxnUKD0F/DVwCLAAuBH4\nkaSDmhrVGEg6DDgZuK/ZsYzRw8AuNY8DmxvOyEmaBdwKCHg7sD/wUWB1M+MapUN5+e/jECCAf2tm\nUKP018BHgI8Bryb9T/4IcPZ4BuFZCWMg6U7gvog4uWbZo8BVETGuv8iiSOoBTouIS5sdy1hJmgGs\nA94ZEVc3O56xkrQGODsiLm52LCMlaXvgv0mJwTnAkog4rblRjZykRcB7IqKULVEVkr4AHBMRv9/s\nWIom6dPAmcCuEfFSs+MZCUnXAM9FxIk1y74D7BQR7xivONxiMEqSOoE3ANfnVl1POmO15ptJ+ow/\n3+xAxiJrXjye1LJzW7PjGaVLSAnzjc0OpAB7Zd1uSyVdIWmvZgc0Cu8E7pR0paTVku6RdJokNTuw\nscji/3Pge2VLCjL/BRwr6dUAkl5Dav386XgGMa79Fi1mNtAOrMotXwW8afzDsTouBO4Bbm92IKMh\n6UBS7FOAHuBdEXF/c6MaOUknA3sD7292LAW4EzgJeAiYA/wNcJukAyLiuWYGNkJ7AacC5wPnAgcD\nX8/WlXL8R2YhsCfwj80OZJS+SDqheUDSFtJ39Ocj4lvjGYQTg7HL98WozjIbZ5LOA44EjoyILc2O\nZ5QeJv3DngX8L+A7kroiYklzwxo+SfuRxuAcFRF9zY5nrCLiZ7WvJd0B/BY4ETivKUGNThuwuKbL\n81eS9iH1Z5c5MTgZ+GVE3NPsQEbpOOADwHuBX5P+/i+UtDQi/mm8gnBiMHrPAluAebnlc9i2FcHG\nkaTzgeOBYyPit82OZ7SyL9LfZC8XSzoU+CtSU2lZHE5qXVtS00rdDhwt6cPA9IjY2KzgxioieiT9\nGtin2bGM0ArggdyyB0mD3UpJ0hzgj0nJTVl9GfhKRFyRvb5f0h6kwYfjlhh4jMEoZf+07yY1XdVa\nSHn7gUtP0oWkbPuNEfFQs+MpWBswudlBjNCPSKP2D655LAauyH4udSuCpCmk0eMrmh3LCN0K7Jdb\nti/weBNiKcpJwEbSZ6usppFOOGttYZy/q91iMDbnAZdJuov0h/ZhYFfgoqZGNULZ6P29s5dtwO6S\nDgbWRMQTzYtsZCR9k9SP/U7geUmV1pyeiOhpXmQjJ+lc4CfAk6Q+x/eSpmOW6loGEbEWWFu7TNJ6\n0merNF0iFZK+AlwNPEFqHfwMMB34TjPjGoXzSWMjPg1cSZpu/THgU02NapSyQYd/AVwRES82O54x\nuBo4S9JSUlfC64EzgO+OZxCerjhGkk4FPkmaP7sE+KuI+EVzoxoZSV3ATXVWfSciThrfaEZP0kAf\n5r+NiEXjGctYSboUOJbUVbWONPf/yxFxXTPjKoKkbso7XfEK4GhS98gzwB3AZyIi3yw/4WUXzPoC\nqeXgCdLYgq9HCb8UJB1LutbH70XEXc2OZ7QkzQT+HngXKfFcQWoB+buI2DBucZTwM2BmZmYN4jEG\nZmZmVuXEwMzMzKqcGJiZmVmVEwMzMzOrcmJgZmZmVU4MzMzMrMqJgZmZmVU5MTAzM7MqJwZmZmZW\n5cTAzBpK0s6SVkg6p2bZQZI2SHpPM2Mzs235kshm1nCS3kK6QcwxwD2kOyzeFREfbGpgZrYNJwZm\nNi4kXQD8EXAzcBRwcNnuemn2u8CJgZmNC0mTgXuBfYAjIuLOJodkZnV4jIGZjZf5wCuBAPZqbihm\nNhC3GJhZw0maBNwOPArcCSwCDoqIJ5oZl5lty4mBmTWcpHOB9wIHAeuAnwFTgWMjYmszYzOzl3NX\ngpk1lKRjgI8DH4iItZHORk4C9gf+upmxmdm23GJgZmZmVW4xMDMzsyonBmZmZlblxMDMzMyqnBiY\nmZlZlRMDMzMzq3JiYGZmZlVODMzMzKzKiYGZmZlV/X8b4vTxVI1W5QAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Lager et plott av den 4. deriverte av 1/(1+16x^2)\n",
"def d4f(x):\n",
" return 6144*(1280*x**4-160*x**2+1)/((1+16*x**2)**5)\n",
"\n",
"x = linspace(0, 8, 1001)\n",
"plot(x, d4f(x))\n",
"xlabel('x')\n",
"title('Fjerde derivert av Runges funksjon');"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"De er dermed ingen overraskelse at feilen er stor og feilestimatet svikter over intervaller som inkluderer $[0,1]$. Dette intervallet må partisjoneres i langt mindre biter for at vi skal få et meningsfyllt resultat. \n",
"\n",
"Men hvordan kan intervallet partisjoneres når $f^{(4)}(x)$ ikke er kjent? "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Adaptiv integrasjon\n",
"Gitt at en grunnleggende integrasjonsrutine som gir en tilnærmelse $Q(a,b)$ med et feilestimat $\\mathcal{E}(a,b)$. \n",
"\n",
"Vi ønsker å finne en partisjon av intervallet $[a,b]$\n",
"\n",
"$$ a = X_0 < X_1 \\cdots < X_m = b $$\n",
"\n",
"slik at\n",
"\n",
"$$ |\\mathcal{E}(X_j, X_{j+1})| \\approx \\frac{X_{k+1}-X_k}{b-a} \\text{Tol} $$\n",
"\n",
"der Tol er en toleranse gitt av brukeren. I så fall vil\n",
"\n",
"$$ \\text{Samlet feil} \\approx \\sum_{j=0}^{m-1} \\mathcal{E}(X_k, X_{k+1}) \n",
" \\leq \\text{Tol}. \n",
"$$\n",
"\n",
"#### Adaptiv integrasjon\n",
"Gitt $f$, $a$, $b$ og Tol.\n",
"* Beregn $Q(a,b)$ og $\\mathcal{E}(a,b)$. \n",
"* **if** $|\\mathcal{E}(a,b)| \\leq \\text{Tol}$:\n",
" * Godkjenn resultatet, returner $Q(a,b) + \\mathcal{E}(a,b)$.\n",
"* **else**:\n",
" * La $c=(a+b)/2$, og prosessen på delintervallene $[a,c]$ og $[c,b]$, med toleranse $\\text{Tol}/2$. \n",
"* De godkjente resultatene fra hvert delintervall summeres underveis. \n",
"\n",
"#### Implementasjon\n",
"Den adaptive algoritmen er implementert med Simpsons metode som basismetode. `simpson_basic` er lett omskrevet fra funksjonen over, den returnerer bare $S2$ og feilestimatet, det er det eneste som behøves i den adaptive algortimen. \n",
"\n",
"`simpson_adaptive` er en rekursiv funksjon (funksjon som kaller seg selv). \n",
"For å unngå at den gjør så i det uendelige, er det lagt til en ekstra variabel `level` som øker med 1 for hver gang funskjonen kalles av seg selv, og en maksverdi for denne. "
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def simpson_basis(f, a, b):\n",
" # Regner ut en tilnærmelse til integralet av f over [a,b].\n",
" # Inn: Funksjonen f, og integrasjonsintervallet [a,b]\n",
" # Ut: S_1(a,b), S_2(a,b) og feilestimatet for S_2(a,b).\n",
" \n",
" # Finner nodene \n",
" c = 0.5*(a+b)\n",
" d = 0.5*(a+c)\n",
" e = 0.5*(c+b)\n",
" \n",
" # Regner ut S1=S_1(a,b), S2=S_2(a,b) og feilestimatet for S2\n",
" H = b-a\n",
" S1 = H*(f(a)+4*f(c)+f(b))/6\n",
" S2 = 0.5*H*(f(a)+4*f(d)+2*f(c)+4*f(e)+f(b))/6\n",
" feilestimat = (S2-S1)/15\n",
" return S2, feilestimat\n",
"\n",
"def simpson_adaptive(f, a, b, tol = 1.e-6, level = 0, maks_level=15):\n",
" # Finner en tilnærmelse til integralet av funksjonen \n",
" # f over intervallet [a,b], gitt toleransen tol.\n",
" \n",
" Q, feil_estimat = simpson_basis(f, a, b) \n",
" \n",
" # -------------------------------------------------\n",
" # Utskrift og plott for demonstrasjonens skyld. \n",
" if level == 0:\n",
" print(' l a b feil_est tol')\n",
" print('==============================================') \n",
" print('{:2d} {:.6f} {:.6f} {:.2e} {:.2e}'.format(\n",
" level, a, b, abs(feil_estimat), tol))\n",
" \n",
" x = linspace(a, b, 101)\n",
" plot(x, f(x), [a, b], [f(a), f(b)], '.r')\n",
" # -------------------------------------------------\n",
" \n",
" if level >= maks_level:\n",
" print('Maksimalt antall nivåer nådd')\n",
" return S2\n",
" \n",
" if abs(feil_estimat) < tol:\n",
" resultat = S2 + feil_estimat # Godkjent, returnerer forbedret approx. \n",
" else:\n",
" # Deler intervallet i to, og bruker metoden på \n",
" # hvert delintervall.\n",
" c = 0.5*(b+a)\n",
" resultat_venstre = simpson_adaptive(f, a, c, tol = 0.5*tol, \n",
" level = level+1)\n",
" resultat_hoyre = simpson_adaptive(f, c, b, tol = 0.5*tol, \n",
" level = level+1)\n",
" resultat = resultat_hoyre + resultat_venstre\n",
" return resultat\n",
" "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"##### Eksempel\n",
"Bruk adaptive Simpson til å finne en tilnærmelse til integralet $ \\int_0^5 1/(1+16x^2)dx $ med ulike toleranser, f.eks. Tol=$10^{-3}, 10^{-5}, 10^{-7}$. Sammenlign den numeriske løsningen med den eksakte. "
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" l a b feil_est tol\n",
"==============================================\n",
" 0 0.000000 8.000000 4.25e-02 1.00e-03\n",
" 1 0.000000 4.000000 1.85e-02 5.00e-04\n",
" 2 0.000000 2.000000 5.11e-03 2.50e-04\n",
" 3 0.000000 1.000000 7.84e-04 1.25e-04\n",
" 4 0.000000 0.500000 6.41e-04 6.25e-05\n",
" 5 0.000000 0.250000 3.43e-05 3.13e-05\n",
" 6 0.000000 0.125000 1.21e-06 1.56e-05\n",
" 6 0.125000 0.250000 1.31e-06 1.56e-05\n",
" 5 0.250000 0.500000 7.82e-07 3.13e-05\n",
" 4 0.500000 1.000000 1.45e-05 6.25e-05\n",
" 3 1.000000 2.000000 1.40e-05 1.25e-04\n",
" 2 2.000000 4.000000 8.29e-06 2.50e-04\n",
" 1 4.000000 8.000000 4.33e-06 5.00e-04\n",
"\n",
"Numerisk løsning = 0.665849, eksakt løsning = 0.384889\n",
"\n",
"Toleranse = 1.0e-03, feil = 2.810e-01\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAe0AAAEACAYAAACebi6nAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3XtgXHWd9/H3d25Jk7Rp6ZUK4S5Q\nQC4tgtyagF14YFdR2a0UFJ5d5YECssuj62VBUHYRn3Wr9cKDsC61ULYouvuIiq1gC0IRaBWh3Aul\nQGlL723umZnv88eZtNNp0kySSc5cPi8cJ3Pmd875fNMk3zlnzjlj7o6IiIgUv0jYAURERCQ/atoi\nIiIlQk1bRESkRKhpi4iIlAg1bRERkRKhpi0iIlIi1LRFRERKhJq2iIhIiVDTFhERKRGxsAPkGjdu\nnB988MEFXWZLSwu1tbUFXWYYyqUOUC3FqlxqKZc6QLUUo6GoY8WKFZvcfXxf44quaR988MEsX768\noMtcunQpjY2NBV1mGMqlDlAtxapcaimXOkC1FKOhqMPM1uQzTrvHRURESoSatoiISIlQ0xYRESkR\natoiIiIlIq+mbWZnmdkvzGytmbmZXZ7HPMeZ2aNm1paZ76tmZoNOLCIiUqHy3dKuA1YC1wFtfQ02\ns1HAb4ENwMnA54AvANcPLKaIiIjk1bTd/dfu/hV3fwBI5zHLJUANcJm7r3T3nwHfBK4fzq1tv/RS\n3IzpTU1gBpdeOlyrFhERKbihek/7Q8Dv3T17q3wRMBk4eIjWuZf0T34OGLteJSxYMFyrFhERKThz\n9/7NYNYMXOPu8/YxZjHwjrv/bda0BmANcJq7P5kz/grgCoCJEydOXbhwYb8y9WZ609kYu+tz4NEl\nSwqy7DA0NzdTV1cXdoyCUC3FqVxqKZc6QLUUo6Goo6mpaYW7T+tr3FBeES331YD1Mh13vxO4E2Da\ntGleqCvNuBnuYDiOYXhJX42nXK4mBKqlWJVLLeVSB6iWYhRmHUO1e3w9MCln2oTM/YYhWufeZl0c\nNG6LkBzbwM6LZg3bqkVERAptqLa0nwS+aWbV7t6emTYDeBd4c4jWuRe7916+0vjXfPKFLsZWTeS1\nES/SNFwrFxERKbB8z9OuM7MTzOyEzDwNmccNmee/YWaPZM1yH9AKzDOzY83s48CXgDne3zfRB+mA\nw45hc0dm435jajhXLSIiUlD57h6fBvwpcxsBfC3z9dczz+8PHNY92N23E2xZTwaWAz8A/g2YU5DU\n/XBiw2jeTW8BoD6633CvXkREpGDy2j3u7kvZfSBZT89f3sO054GzBhqsUI4/YDT3VxmNQH1iP5LJ\nJLFY0X0iqYiISJ/K/trjtVUxnh9zNB2pNqqiI1i99KmwI4mIiAxI2TdtgIPG17GjaxsA7zzxp5DT\niIiIDExFNO3DRkfY2bUdgNRmHYwmIiKlqSKa9uGjo+xMBVvaVenakNOIiIgMTEU07Yk1xrbITgBq\nYyNDTiMiIjIwFdG0zYxtoxMA1Klpi4hIiaqIpg2QOHEaaU9RExvJjo0bw44jIiLSbxXTtI8+/iRa\nks0AvPbQoyGnERER6b+KadrHN4yhNdO0t776dshpRERE+q9imvbI6jitqaBpp7aHHEZERGQAKqZp\nA7Smg6adSFeHnERERKT/Kqppd8Y6ABgR0bnaIiJSeiqqadvoKAA1MTVtEREpPRXVtOuPPACAEWra\nIiJSgiqqaU+ZcTZpTzMiWsuO9evDjiMiItIvFdW0x06YQFvmCPIXFz8WchoREZH+qaimDdCWbAVg\n62trQ04iIiLSP5XXtFNB007vSIacREREpH8qrmm3p4OmHU9WhZxERESkfyqvaXvQtKvRBVZERKS0\nVFzT7ooHF1ipjtaEnERERKR/Kq5pe31wgZURatoiIlJiKq5pjzx0MgDVMTVtEREpLRXXtA+bfiZp\nT1MdqaFly9aw44iIiOSt4pr2+xoa6Ei1Yma8/vsnwo4jIiKSt4pr2rFoZNe52hteeCPkNCIiIvmr\nuKYN0J5qA6Bzc3vISURERPJX0U071hELOYmIiEj+KrJpd3jQtOOuC6yIiEjpyLtpm9lsM1ttZu1m\ntsLMzuxj/Cwze9bMWs1svZnda2aTBh958DoiwW7xahsRchIREZH85dW0zWwmMBe4FTgRWAY8ZGYN\nvYw/HbgH+DFwDHAhMAVYUIDMg5asSgNQFVXTFhGR0pHvlvb1wDx3v8vdX3L3a4F1wFW9jP8Q8I67\nf9vdV7v7H4DvAacMPvLgjd+5lklzZnDSdy+GqipYujTsSCIiIn3qs2mbWQKYCizOeWoxcFovsz0B\n7G9mf2WBccAngV8PJmyhnDLvh0SSnRhAZyece27YkURERPpk7r7vAWaTgbXAdHd/LGv6V4FL3P3I\nXub7BHA3MAKIAb8FPuqeOQpsz7FXAFcATJw4cerChQsHVk0vmpubqaur2/V4elNT0LAzHHh0yZKC\nrnMo5NZRylRLcSqXWsqlDlAtxWgo6mhqalrh7tP6HOju+7wBkwn62pk5028CXu5lnikEjf4LwAeA\nc4HngPl9rW/q1KleaEuWLNnjcTqR8JSZO3ga3BOJgq9zKOTWUcpUS3Eql1rKpQ531VKMhqIOYLn3\n0R/dPa/3tDcBKSD3yO8JwIZe5vky8LS7/6u7P+fui4DZwKfM7MA81jmkbNEittWPJm1GR1UVLFoU\ndiQREZE+9Xl1EXfvNLMVwAzgp1lPzQB+1stsNQSNPlv3YyNsjY3837/9DF0jaxjXvJNrGhvDTiQi\nItKnfC8JNge4x8yeJjjI7EqC3eZ3AJjZfAB3/3Rm/IPAXWZ2FbAI2B/4DvBHd3+rcPEHLpIKTvvq\n6EyHnERERCQ/eTVtd7/fzMYCNxA04JXA+e6+JjOkIWf8PDMbCVwD/BuwHVgC/GOhgg9WNB1s+Cct\n/A1/ERGRfOR98W13vx24vZfnGnuY9j2Cc7OLUsyCo+bTkYq8kquIiJSgiu1YVVXB65VULBpyEhER\nkfxUbNMeOW4MAOloxX4LRESkxFRsx2o4agoAri1tEREpERXbtE+afg6447EYW9a/G3YcERGRPlVs\n0x5VPwpLJgH409LfhZxGRESkbxXbtAEiyeC0r3dWrQo5iYiISN8qu2lnLrCyc/vOkJOIiIj0rbKb\ndmZLu7Mz94qrIiIixaeim3Y0HWxpJyv72yAiIiWiortVLHMF05SuiiYiIiWgorvVrqui6QIrIiJS\nAiq6W9WPrQcgrQusiIhICajopn3QlKMAXRVNRERKQ0U37RPPmrHrqmjvrVsbdhwREZF9quimXVdf\nv/uqaEseCTmNiIjIvlV004bdF1hZ+/rrIScRERHZNzXtzAVWWrbpqmgiIlLc1LRTQdPuSKZDTiIi\nIrJvFd+0Y5nd4ymv+G+FiIgUuYrvVLHMdyCtC6yIiEiRq/hONaI6AUAqqnO1RUSkuFV80x47aSwA\n6VjFfytERKTIVXynOuL4EwFIx6N06mA0EREpYhXftI859UzwNERjvPLi82HHERER6VXFN+34iGos\nc672y8seCzmNiIhI7yq+aQNEuoKmvfmdd0NOIiIi0js1bXZfYKWtpT3kJCIiIr1T0waimQPQMp8d\nIiIiUpTybtpmNtvMVptZu5mtMLMz+xifMLOvZ+bpMLO3zOxzg49ceDEyV0WL6DWMiIgUr1g+g8xs\nJjAXmA08nrl/yMymuPtbvcz2n8CBwBXAa8BEYMSgEw+BRMwAXRVNRESKW15NG7gemOfud2UeX2tm\n5wFXAV/OHWxmfwF8GDjM3TdlJr85yKxDpn5ULVscUjFdFU1ERIpXn5uWZpYApgKLc55aDJzWy2wX\nAs8A15vZO2b2mpl918zqBpV2iBx4xKFAcIGVVNpDTiMiItIzc993kzKzycBaYLq7P5Y1/avAJe5+\nZA/z/AZoBB4Bvg6MBr4HPOfuF/Uw/gqC3ehMnDhx6sKFCwdaT4+am5upq+v99ULrtk08/exKcOfY\nkz7IuPragq6/UPqqo5SoluJULrWUSx2gWorRUNTR1NS0wt2n9TUu393jALnd3XqY1i2SeW6Wu28H\nMLNrgEVmNtHdN+yxYPc7gTsBpk2b5o2Njf2I1belS5fS1zKfWf4sHotR3baVxo9eUND1F0o+dZQK\n1VKcyqWWcqkDVEsxCrOOfI682gSkgEk50ycAG/YeDsA6YG13w854KXPf0K+Ew8Qy53utf/XVkJOI\niIj0rM+m7e6dwApgRs5TM4Blvcz2BDA55z3s92fu1/Q35HDoPle7edvOkJOIiIj0LN9znOYAl5vZ\nZ8zsaDObC0wG7gAws/lmNj9r/H3AZuBuMzvGzE4nOGXsAXd/r4D5Cyaauf54MqkD0UREpDjl9Z62\nu99vZmOBG4D9gZXA+e7evdXckDO+2cw+THDw2TPAVuC/gS8VKnihxdKZC6yYztUWEZHilPeBaO5+\nO3B7L8819jDtFeAvBpxsmHVfYCUV1bnaIiJSnLRZmVFXH7z9rgusiIhIsVLTzjj4yCOA4AIrfZ27\nLiIiEgY17YypHz4P3PFYjLVvvxN2HBERkb2oaWfUjxsfnKttxp+W/CbsOCIiIntR084SyZz2tW7V\nGyEnERER2ZuadpZoV9C0W3a2hZxERERkb2raWaKpoGl3pS3kJCIiIntT084Sy3z+SSqib4uIiBQf\ndacs1VXBtWZSsf58+JmIiMjwUNPOMm7yBCA4V1tERKTYqGlnmXb2OQCk4zFa21pDTiMiIrInNe0s\nh618mRtuuYWbbrmFETW1cN99YUcSERHZRU072yWXEPU0lvVYRESkWKhp57DMZcd10peIiBQbNe0c\nnnMvIiJSLNS0sy1YQNfICB6BjrFRWLAg7EQiIiK76ITkbLNmsWj7MmqOfIi2dQfyl7NmhZ1IRERk\nF21p59j/kBMAiI3cSWcyHXIaERGR3dS0cxx96sdwN2I123n19VfCjiMiIrKLmnaO2tFjSbaMxiLO\n60//Iuw4IiIiu6hp9yC5cyQAbRtfDTmJiIjIbmraPUjuHAGA+baQk4iIiOympt2DdHvQtOPV7SEn\nERER2U1Nuwex2GgA4iPVtEVEpHioafdgwvumABAftZP2rlTIaURERAJq2j04+vSLgtO+6rby6uo3\nw44jIiICqGn3qH7CQSRb6rFImlVPPhB2HBEREUBNu1fJ7fUAdG56MeQkIiIiATXtXnRlTvuKsCPk\nJCIiIoG8m7aZzTaz1WbWbmYrzOzMPOc7w8ySZrZy4DGHX7qzDoBErY4gFxGR4pBX0zazmcBc4Fbg\nRGAZ8JCZNfQx3xhgPvDIIHMOuxG1kwGI17eGnERERCSQ75b29cA8d7/L3V9y92uBdcBVfcz3I+DH\nwJODyBiKQ4//MADx0VvY3qzGLSIi4euzaZtZApgKLM55ajFw2j7mmw1MAv55MAHDcvhJ55LqrCZa\n1crzf3gw7DgiIiKYu+97gNlkYC0w3d0fy5r+VeASdz+yh3mOAx4GTnX31WZ2M3CRux/byzquAK4A\nmDhx4tSFCxcOsJyeNTc3U1dX1+/5WjbcSPX4d1m7/EwaPnh5QTMNxEDrKEaqpTiVSy3lUgeolmI0\nFHU0NTWtcPdpfY2L9WOZud3depiGmVUBC4HPu/vqvBbsfidwJ8C0adO8sbGxH7H6tnTpUgayzP/3\nwzoYD9XR5gHNX2gDraMYqZbiVC61lEsdoFqKUZh15NO0NwEpgl3d2SYAG3oYvz8wBbjbzO7OTIsA\nZmZJ4Hx3z93VXpS6WmsBiNe0hZxEREQkj/e03b0TWAHMyHlqBsFR5LnWAscBJ2Td7gBWZb7uaZ6i\nFI1PACAxujnkJCIiIvnvHp8D3GNmTwNPAFcCkwmaMWY2H8DdP+3uXcAe52Sb2XtAh7uX1LnaBx17\nDlvTi4iP3kJnRyeJqkTYkUREpILldcqXu98P/D1wA/AscAbBbu41mSENmVtZ+cDpHyHVUUO0qpWn\nlxT24DgREZH+yvuKaO5+u7sf7O5V7j41+0hyd29098Z9zHtzb0eOF7NYPE7nlrEArHv50ZDTiIhI\npdO1x/vQtS04rN98c8hJRESk0qlp98G7gk/7SozSVdFERCRcatp9GPu+4wGI77edzmQ65DQiIlLJ\n1LT7MG3G3+GpKPFRm1nxjN7XFhGR8Khp96F29Fg6t43HzFnz9M/CjiMiIhVMTTsPnVuC97UjybUh\nJxERkUqmpp0HbxsFQJU+W1tEREKkpp2H8fsHp5gnxm2luSMZchoREalUatp5OOm8/0U6GScxajNP\nPfGbsOOIiEiFUtPOQ039eDq3BB8e8t6f/ivkNCIiUqnUtPPUuTk4GC0e3RhyEhERqVRq2nmydLCl\nXT1uZ8hJRESkUqlp52nKqZ8AIDF+A++u06lfIiIy/NS08/T+k8+na+d+ROMdPP2r74cdR0REKpCa\ndj90bAg+prNr+6shJxERkUqkpt0PXS1jAKjab0fISUREpBKpaffD/kecDUD1pPW0tTSHnEZERCqN\nmnY/nHLu5XQ1jyZa1coj938z7DgiIlJh1LT7IRaP07E+OPWrbdvKkNOIiEilUdPup1Rb0LRHjN+O\nu4ecRkREKomadj8d33g5AFWT1vL8cyvCDSMiIhVFTbufDju+iY5N+xOJJnnh4R+EHUdERCqImvYA\ntK8fB0B1zYaQk4iISCVR0x6A/caeAMCIA9azYatO/RIRkeGhpj0Ap3/iiyRbRxGv3c6jP/mXsOOI\niEiFUNMegHj1CNrengxAtOOFkNOIiEilUNMeoETs/QDUHLiB5tb2kNOIiEglUNMeoLMvvplUey2J\n+k08slC7yEVEZOjl3bTNbLaZrTazdjNbYWZn7mPsx81ssZltNLOdZvaUmX2kMJGLQ3VdPa1vHQhA\nqvnZkNOIiEglyKtpm9lMYC5wK3AisAx4yMwaepllOvA74ILM+F8D/7WvRl+KYtGjAag5ZK0+QERE\nRIZcvlva1wPz3P0ud3/J3a8F1gFX9TTY3a9z99vc/Wl3X+XuXwNWABcWJnZx+ItPfZ2u5jHEa7fz\n6//4QthxRESkzPXZtM0sAUwFFuc8tRg4rR/rGgls7cf4opeorqHmoRhTP/YeH7vuh2AG990XdiwR\nESlT+WxpjwOiQO7lvzYAk/JZiZldDRwA3NOvdCXg1B+tYOTOFiLuOMAll4QdSUREypT19UlVZjYZ\nWAuc5e6/z5p+E3Cxux/Vx/yfIGjWn3T3X/Qy5grgCoCJEydOXbhwYb+K6EtzczN1dXUFXWa36U1N\nWNZjBx5dsmRI1jWUdQw31VKcyqWWcqkDVEsxGoo6mpqaVrj7tD4Huvs+b0ACSAJ/nTP9B8Cjfcz7\nCaAVuKiv9XTfpk6d6oW2ZMmSgi9zF/AUuIMnzTxtNmSrGtI6hplqKU7lUku51OGuWorRUNQBLPc8\nemSfu8fdvZPgILIZOU/NIDiKvEdm9jfAvcDl7v5An68eStWCBZgZqUiEVQ2H8P2rLgs7kYiIlKl8\njx6fA1xuZp8xs6PNbC4wGbgDwMzmm9n87sFm9klgAfAl4DEzm5S57Vfg/OGbNQtLp7l+7teZfvfP\nuO9Df0lru66QJiIihZdX03b3+4G/B24AngXOAM539zWZIQ2ZW7crgRjwHYJTw7pvPy9M7OJz48V/\nR/3O7ax+32F8a963wo4jIiJlKO8rorn77e5+sLtXuftUd38s67lGd2/MeWw93Bp7WnY5GDd2Emet\nfhyAxTXvJ9nVFXIiEREpN7r2eAHdeNHFjGrewaoD388td30j7DgiIlJm1LQLqOF9h3LO6mAHxC/3\nO4n1mzaGnEhERMqJmnaB3XbZ1UzcvIG1Ew/g5vnfDzuOiIiUETXtAquvH8PfbPszAL89cga/efy3\nIScSEZFyoaY9BL50+fUc+8afaamp4/vPryKdSocdSUREyoCa9hCIxmLccPREqjraWX7Uh7jh9q+H\nHUlERMqAmvYQaTztPP7q1WDX+E8PPocly3I/JE1ERKR/1LSH0Jwr/5Fj3niOnXUj+dorG9m0dXPY\nkUREpISpaQ+hRKKKb558OKN3bOPlg4/hmvn/Tjq9709VExER6Y2a9hCb9oHTuK5tJdFUkqUfOJfr\n5t4UdiQRESlRatrD4KpPXsPMlx8C4GfHfZSbbr8l5EQiIlKK1LSHyZxrbmTG84tIR6Pcfdi53HrX\nbWFHEhGREqOmPYzuvup6znhhCZ2JBHc0nM1Xb//nsCOJiEgJUdMeRrF4nPuumM3pLyylM5HgR+8/\nj+u+cxPuOjhNRET6pqY9zBKJKn5y5TV8+PlFpKIx7j/+Y1z63W+yo3l72NFERKTIqWmHIBqLce/n\nvsjMlQ8SSaV45APn8ZGf/pLfPvGbsKOJiEgRU9MO0dxrb+QL7z626zzu2Vtr+cfv3kyyqyvsaCIi\nUoTUtEP2D5/+B+49JB5cOa12JPOPu5Dz593Dg4/8POxoIiJSZNS0i8C0D5zO4stmMXPlg9S0tfLc\n4SdxdfIALv3ON3j+5efCjiciIkVCTbtIRGMx5l57I/PHd3Dyy8voTCR4+Pj/wcdW7eBbV1yGRyJM\nb2oCM7jvvrDjiohICNS0i8wZJzfx4FWzueW9ZRy55kWaa+u4dt5/4u4Y4ACXXBJyShERCYOadpH6\n7MzZ/O5TM/n8m4+Q6Ora9Q9lQMqMq3/4Td5+b12YEUVEZJipaRexaDTK5//n/969hU3QsF/f7wB+\ntfpYzvr2k5x6y218/u45vPPeu2FGFRGRYRALO4DkYcEC7JJLgsZtxoK/Po+aqrdp7TiQ9S3H8cAr\n8MCrTzGmejVH1Tdz4YkncdHp5xGN6Z9XRKSc6K96KZg1C2bN4tGlS2lsbORrwNeAXz31CD96Yhkv\nbx9Pa8eBbG07kifb4MmH4CuLf8q46rc4bGQn0484nJnTL2B03aiwKxERkUFQ0y5hF5xyDheccg4A\nf3jxj/zo0UU8tynKxvYGUqlRbGg5lg0tsGw9fOPx31GTWMeEqs0cNMqYeuABnH9yI4dPPijkKkRE\nJF9q2mXi1CknceqUkwBIJZM88MQifv3cn3lte4SNHZPo6hpPa0cDb3Y08OYOePQdmPPkSmKx31MX\n38R+iWYm1TiHjR3DiQcfQdPxpzBm5OiQqxIRkWxq2mUoGosxc/oFzJx+wa5pL721ip8/+TDPrt3I\nO60JtnaOpb1rPMnkGLYlx7CtDd7YDsvWwT0rk/DL3xON7WBEdBu1sWbqEx2MrYqw/6gRHLTfeI46\n4GCmHn4M40aPDbFSEZHKoqZdIY5uOJx/ajh8j2ktbS0s+uPj/OG1l3hjy042tEXZ1jWStuQYksl6\nUsnRNCdH09wBG1oyM63vnnsn8Acs0ko82kxVtIXqSBsjop3UxVOMjBv11XHG1oxg/MhRTN5vLAeO\nm8yUZcsZ/b+uYnr3YhYsCN6zFxGRPuXdtM1sNvAFYH/gBeDv3f33+xg/HZgDHAO8C/wfd79jcHGl\nkGpH1PLx08/l46efu9dzO1p28vgLz/Dsm6t4c/MW3mvpZGtnlJ3JEbSlaulMjiSVqsPTNXSma+js\nCtr4vm3glX+9FgzMwQ24+lJaX/s8yUiCZKSKZKSaVHQEqXgN6Vgt6XgtXlUHVaOgup5I9SiitWOI\n1IwhXlNPrGYMiZH7UT1yLLGqmiH4LomIFI+8mraZzQTmArOBxzP3D5nZFHd/q4fxhwC/Bv4DuBQ4\nA7jdzDa6+88KFV6GzqjakZz/wbM5/4Nn9zqmq6uTF99exYtrVrFm43re3bGdLa0d7OhI05yM0JaK\n055K0Jmupis9gmS6mkQ6iWXmNwe2O7XeAqkWSA0us2dOZnczHCON4URIEyFt0d33Ftt9H4ll7uO4\nxfFInHQkhkfjeCSBRxN4JA6xBEQSEI1DrAqiVcG0aILYhk2saX8ei1dhsargPlpFJF5FJFaNxRNE\n49VYLEEsXgWxKqLxBNF4FdFYFbHECCwaHVzxg3XffXDJJdoDIlLk8t3Svh6Y5+53ZR5fa2bnAVcB\nX+5h/JXAu+5+bebxS2Z2CvB5QE27TMTjCY4/dArHHzol73n81otx2H3BmCi8edKX8JYt0L4V69hB\ntHMnka4Woql2oqk2oqlOot5JNN1FlCRRTxEhRQTHcCzTrS3zasBwwAnaYOaVQPfVaZyCOxRgkBen\n8x5yuQGZlziZCvdxn/u1ge09ffc0w80yyzdit62hxiDikDZ475rL+HTqMcwiwQiLYFjWvRHJXJup\ne1rEgseRrDHBdAMiRDKPdy/TssZFMCPrPncdtscyzXY/nz2++78NG9azePHrwV6dnCyW+UHpnr47\nM3tkAYjsGtO9DPZY/+4xwK71B19HMt/vXVczzCw7suvn1HZNz15+xLKnw6qNq1j7Z7Kes5x5gwwA\nFrE9lh2xnh93/1xFehu/a9zu78Me68sZ1z3/rvG5682s58Ud64mveamHfHte52t3rj1l/3vl2itz\n7veZPWvf9ffC9lzX7r8jOZnMSCxbxvtmzmJ6eyfpqgSR3yyCxsa9sgylPpu2mSWAqcC3cp5aDJzW\ny2wfyjyfbRFwmZnF3V0fGF2hbMECyFwoxgB+vIBDPjL4LTpPpejsaKFj52a6WreRbN1Oqn0nqbYd\npDtaSLc3412teGcrdLVCsh2SHViyHVKdRFIdWLqTSKoLSyexdBfmSSKexDxFJJ3EPE2EVPDY08G2\nvKeDx7tb4R43gla5Z5vNadBmu/9Q7DE9qCy7yn5+U/KcBqTb0kQyz0Ucxm1Psi79RP/WV2RWlNNV\nfp8NO0ABLQ07wMD911dewzo6g42Ojk4491zo6BjWDPlsaY8DosCGnOkbgA/3Ms8k4OEexscyy9vj\n18nMrgCuAJg4cSJLly7NI1b+mpubC77MMJRFHZMnw5IlNDc3U1dXF0wbsppGBrcoUDtEq4A9axkA\nT6VwkngqCakkpJOQDl48kE7t+to9BekUlk6BpyAdvKDovrmnsHR612PS3S8sUuBp8BTmvsc0cwdP\ncXzsPuJdTgRIA61VET4YuRh3p/u/NJ71GNKkIfuxBy9SgpHsMW/wTDrzmiEzzXd/Tdb/9zSdnFFZ\nSyH43+5dKQ6k06ldW509zZukevqHAAAHdklEQVQ9fe97MuP2nkp3Psse6btnMd/r8R7juh9a7rL3\nfDW1RxL3zNZgT3n2njffcX0vZ8+sfc+/7zx9j8tZruUM6yVHf5bd97z7/h4csr6DaPeLW8A7O3l0\nuP8mu/s+b8BkgkrOzJl+E/ByL/O8CtyYM216ZjmT9rW+qVOneqEtWbKk4MsMQ7nU4a5ais6CBe6Z\n1u4QPC5hZfFvkqFaikgiEfx+dN8SiYItGljuffRjd8/rA0M2EbwxOCln+gT23vrutr6X8Ulgcx7r\nFJHhNGsWuPPokiXBnyMdhCayt0WLIJEItscTieDxMOuzabt7J7ACmJHz1AxgWS+zPcneu85nELyS\n0PvZIiJSehoboaMjeHHb0THsB6FB/h/NOQe43Mw+Y2ZHm9lcgt3mdwCY2Xwzm581/g7gADP7Tmb8\nZ4DL2ftgNhEREclTXqd8ufv9ZjYWuIHg4iorgfPdfU1mSEPO+NVmdj7wbYLTwt4FPuc6R1tERGTA\n8r4imrvfDtzey3ONPUx7FDhpwMlERERkD/nuHhcREZGQqWmLiIiUCDVtERGREmHe00WPQ2RmG4E1\nfQ7sn3EE55uXunKpA1RLsSqXWsqlDlAtxWgo6jjI3cf3NajomvZQMLPl7j4t7ByDVS51gGopVuVS\nS7nUAaqlGIVZh3aPi4iIlAg1bRERkRJRKU37zrADFEi51AGqpViVSy3lUgeolmIUWh0V8Z62iIhI\nOaiULW0REZGSp6YtIiJSIsq6aZvZbDNbbWbtZrbCzM4MO1N/mdlZZvYLM1trZm5ml4edaaDM7Mtm\n9oyZ7TCzjWb2oJkdG3au/jKzq83suUwdO8zsSTO7IOxchWBmX8n8nH0/7Cz9ZWY3Z7Jn39aHnWug\nzGx/M/tx5nel3cxeNLPpYefqDzN7s4d/EzezX4Wdrb/MLGpmt2T1lNVm9s9mlvdneBRC2TZtM5sJ\nzAVuBU4k+Ozvh8ysYZ8zFp86gk9Vuw5oCznLYDUSfOjMacDZQBJ42Mz2CzPUALwDfJHgA3GmAb8D\n/tvMPhBqqkEys1OBzwLPhZ1lEF4h+CTC7ttx4cYZGDMbDTwBGHABcDRwLfBemLkG4GT2/Pc4CXDg\nJ2GGGqAvAlcDnwOOIvibfDXw5eEMUbYHopnZU8Bz7v7ZrGmvAQ+4+7B+kwvFzJqBa9x9XthZCsHM\n6oDtwIXu/mDYeQbDzLYAX3b3H4adZSDMrB74I0HT/iqw0t2vCTdV/5jZzcBF7l5ye29ymdmtwHR3\nPz3sLIVkZv8EfAGY7O6tYefpDzP7JbDZ3S/LmvZjYKy7/+Vw5SjLLW0zSwBTgcU5Ty0m2MqT4jCS\n4Gdwa9hBBiqzy+yTBHtEloWdZxDuJHhB+7uwgwzSoZm3klab2UIzOzTsQAN0IfCUmd1vZu+Z2bNm\ndo2ZWdjBBiqT/e+Ae0utYWc8DjSZ2VEAZjaFYI/hr4czxLDuix9G44AosCFn+gbgw8MfR3oxF3gW\neDLsIP1lZscR5K4GmoGPufvz4aYaGDP7LHA48KmwswzSU8DlwMvABOAGYJmZHePum8MMNgCHArOB\nbwO3AScA38s8V3LHG2TMAA4B/j3sIAP0TYINjRfNLEXQP//F3W8fzhDl2rS75e77tx6mSQjMbA5w\nBnCGu6fCzjMArxD8IR0NfAL4sZk1uvvKcGP1j5kdSXDcx5nu3hl2nsFw94eyH5vZH4A3gMuAOaGE\nGrgIsDzrrbw/mdkRBO+hlmrT/izwjLs/G3aQAZoJfBqYBbxA8Ps/18xWu/uPhitEuTbtTUAKmJQz\nfQJ7b33LMDOzbwOfBJrc/Y2w8wxEpsGtyjxcbmYnA/9AsPuvlHyIYM/Uyqw9r1HgLDO7Eqh1946w\nwg2Guzeb2QvAEWFnGYB1wIs5014iOPip5JjZBOCjBC86StW/At9y94WZx8+b2UEEB6INW9Muy/e0\nM39QVxDsjsk2g9J+37HkmdlcgleqZ7v7y2HnKaAIUBV2iAH4b4IjrE/Iui0HFma+LtmtbzOrJjjK\nd13YWQbgCeDInGnvp/AfWzxcLgc6CH6uSlUNwcZgthTD3EfLdUsbgt1h95jZ0wS/AFcCk4E7Qk3V\nT5kjrA/PPIwADWZ2ArDF3d8KL1n/mdkPCN43vRDYambde0Ka3b05vGT9Y2a3Ab8C3iZ4j2sWwels\nJXeutrtvA7ZlTzOzFoKfr1Lb1f8t4EHgLYK9ajcCtcCPw8w1QN8meD/+n4D7CU5b/RzwlVBTDUDm\nALTPAAvdfWfYeQbhQeBLZraaYPf4icD1wPxhTeHuZXsjOJDjTYJXeCuAs8LONIAaGgneh8+9zQs7\n2wBq6akOB24OO1s/65hHsMXTQXDe7MPAuWHnKmB9S4Hvh51jALkXAu8S7B1YC/wMmBJ2rkHUcwHw\nZ6AdeJWgaVvYuQZQR1Pm9/yDYWcZZB0jge9kfvfbCI6XuBWoHs4cZXuetoiISLkpy/e0RUREypGa\ntoiISIlQ0xYRESkRatoiIiIlQk1bRESkRKhpi4iIlAg1bRERkRKhpi0iIlIi1LRFRERKxP8HR5h1\nNKve6zAAAAAASUVORK5CYII=\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"def runge(x):\n",
" return 1/(1+(4*x)**2)\n",
"\n",
"a, b = 0, 8\n",
"I_eksakt = 0.25*(arctan(4*b)-arctan(4*a))\n",
"\n",
"tol = 1.e-3\n",
"resultat = simpson_adaptive(runge, a, b, tol=tol)\n",
"\n",
"# Utskrift\n",
"print('\\nNumerisk løsning = {:8f}, eksakt løsning = {:8f}'\n",
" .format(resultat, I_eksakt))\n",
"feil = I_eksakt - resultat\n",
"print('\\nToleranse = {:.1e}, feil = {:.3e}'.format(tol, abs(feil)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Andre kvadraturformler\n",
"\n",
"Simpsons metode er bare et eksempel på kvadraturformler, og alt som er gjennomgått for Simpsons metode kan også gjøres for andre. \n",
"\n",
"\n",
"\n",
" \n",
"La meg til sist nevne noen vanlige klasser av metoder:\n",
"\n",
"**Newton-Cotes metoder** Basert på jevnt fordelte noder over intervallet, se \n",
"[her for noen eksempler](https://en.wikipedia.org/wiki/Newton–Cotes_formulas). Trapes og Simpson er eksempler på lukkede Newton-Cotes metoder, midtpunkt-metoden et eksemp på en åpen Newton-Cotes metode. \n",
"\n",
"**Gauss-Legendre kvadratur**: Velg nodene som nullpunktene i polynomet\n",
"$$ L_m(t) = \\frac{d^m}{dt^m}(t^2-1)^m. $$ \n",
"Disse har presisjonsgrad $d=2m-1$, og er det beste som kan oppnås med $m$ noder. Midtpunkt-metoden er et eksempel på et Gauss-Legendre kvadratur med $m=1$. \n",
"\n",
" \n",
"\n",
" \n",
"\n",
"\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Appendix\n",
"\n",
"#### Den generaliserte middelverdisetningen\n",
"\n",
"_La $f\\in C[a,b]$. Gitt $k$ noder $x_i$ i $[a,b]$ og $k$ positive vekter $w_i>0$, $i=1,\\dotsc,k$. Da fins en \n",
"$\\eta \\in [\\min{x_i},\\max{x_i}]$ slik at_\n",
"\n",
"$$ \\sum_{i=1}^{n} w_i f(x_i) = f(\\eta) \\sum_{i=1}^k w_i. $$"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
}
],
"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.6.3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}