{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "568bde50",
   "metadata": {},
   "source": [
    "# Innlevering 4 - Numerisk integrasjon"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2218e4e2",
   "metadata": {},
   "source": [
    "Standardnormalfordelingen er gitt ved \n",
    "$$\n",
    "    f(x) = \\frac{1}{\\sqrt{2\\pi}} e^{-x ^2/2}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ef54163",
   "metadata": {},
   "source": [
    "Standardnormalfordelingen er et eksempel på en sannsynlighetstetthetsfunksjon, \n",
    "og for TMA4400 holder det å vite at sannsynligheten for å havne mellom $a$ og $b$ er gitt ved \n",
    "$$\n",
    "\\frac{1}{\\sqrt{2\\pi}} \\int_a^b e^{-x ^2/2} \\; dx.\n",
    "$$\n",
    "Det er bare ett problem - standardnormalfordelingen kan ikke antideriveres. \n",
    "Under er kode for å beregne arealet ved riemannsummer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25ccb8a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "#integrasjonsgrenser\n",
    "a=0\n",
    "b=1\n",
    "\n",
    "#gitter\n",
    "n=1000\n",
    "x=np.linspace(a,b,n+1)\n",
    "h=b/n \n",
    "\n",
    "#venstre riemannsum\n",
    "f=h*(np.sum(np.exp(-x**2/2))-np.exp(-b**2/2))/np.sqrt(2*np.pi)\n",
    "\n",
    "print(f)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b000b17a",
   "metadata": {},
   "source": [
    "1. Modifiser koden til å beregne arealet ved trapesmetoden, \n",
    "kjør alt med $a=0$, $b=1$, \n",
    "sammenlikne begge metoder med den korrekte sannsynligheten\n",
    "$$\n",
    "\\frac{1}{\\sqrt{2\\pi}} \\int_0^1 e^{-x ^2/2} \\; dx \\approx 0.3413447460685429485852.\n",
    "$$\n",
    "Kan du forklare forskjellen? (Fun fact: Det er denne beregningen som forteller deg at om du har IQ på 115 er du smartere enn om lag 84 prosent av befolkningen. Mer om dette i TMA4240/4245.)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b05cee1",
   "metadata": {},
   "source": [
    "2. Man kan bruke polynominterpolasjon til å finne bedre teknikker for numerisk integrasjon. \n",
    "Utled Simpsons metode ved å finne lagrangepolynomene til gitteret\n",
    "$$\n",
    "    \\left\\{ a, \\frac{a+b}{2}, b \\right\\}\n",
    "$$\n",
    "og integrere dem.\n",
    "Test Simpsons metode på arealet over og forklar forskjellen.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ccb45b23",
   "metadata": {},
   "source": [
    "De numeriske integrasjonsrutinene gir oss et sett med korresponderende teknikker for å løse ordinære differensiallikninger numerisk. \n",
    "La oss si vi ønsker å løse initialverdiproblemet\n",
    "$$\n",
    "    \\dot x = f\\left( x \\right) \\hspace{10mm} x(0) = 0\n",
    "$$\n",
    "på intervallet $[0,T]$.\n",
    "Vi deler tidsintervallet inn i en partisjon med $n+1$ punkter\n",
    "$$\n",
    "    t_k = kh \\hspace{10mm} 0\\leq k \\leq n\n",
    "$$\n",
    "der $h=T/n$,\n",
    "og integrerer differensiallikningen\n",
    "fra ett gitterpunkt til neste:\n",
    "$$\n",
    "    x\\left( t_{k+1} \\right) - x\\left( t_k \\right) = \\int_{t_k}^{t_{k+1}} f\\left( x\\left( t \\right) \\right) \\; dt\n",
    "$$\n",
    "Hvis vi setter inn en enkelt rektangulær riemannboks (venstre, høyre eller sentrert) eller trapesmetoden for integralet på høyresiden,\n",
    "og så setter inn approksimasjonen $x_k \\approx x\\left( t_k \\right)$ overalt,\n",
    "får vi nå fire metoder på rappen:\n",
    "\\begin{align*}\n",
    "    \\text{Fremlengs Euler:} \\hspace{10mm} &x_{k+1} = x_k + hf\\left( x_k \\right) \\\\[3mm]\n",
    "    \\text{Baklengs Euler:} \\hspace{10mm} &x_{k+1} = x_k + hf\\left( x_{k+1} \\right) \\\\[3mm]\n",
    "    \\text{Midtpunktmetoden:} \\hspace{10mm} &x_{k+1} = x_k + hf\\left( \\frac{x_k + x_{k+1}}{2} \\right) \\\\[3mm]\n",
    "    \\text{Trapesmetoden:} \\hspace{10mm} &x_{k+1} = x_k + \\frac{h}{2}\\left( f\\left( x_k \\right) + f\\left( x_{k+1} \\right) \\right)\n",
    "\\end{align*}\n",
    "Under finner du kode for fremlengs Euler\n",
    "på det logistiske problemet\n",
    "$$\n",
    "\\dot{x}=a\\left( 1-x/k \\right)x \\hspace{10mm} x(0) = x_0 \\hspace{10mm} a >0,\\; k > 0.\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2048adf",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "#parametre\n",
    "a=1 #vekstfaktor\n",
    "k=1 #bærekapasitet\n",
    "\n",
    "#tidsakse\n",
    "T=10\n",
    "N=10000\n",
    "h=T/N\n",
    "t=np.linspace(0,T,N+1)\n",
    "\n",
    "#løsning\n",
    "x=np.zeros(N+1)\n",
    "x[0]=.1\n",
    "for i in range(N):\n",
    "    x[i+1]=x[i]+h*x[i]*(1-x[i])\n",
    "\n",
    "plt.plot(t,x)    \n",
    "plt.savefig('logistisk.png') "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "85686636",
   "metadata": {},
   "source": [
    "3. Implementér de forskjellige metodene for det logistiske problemet \n",
    "og sammenlikne med den analytiske løsningen\n",
    "$$\n",
    "x(t)=\\frac{kx_0e^{at}}{k-x_0+x_0e^{at}}.\n",
    "$$\n",
    "Kan du forklare forskjellene?\n",
    "(Velg vekstfaktoren $a$ og bærekapasiteten $k$ selv. For midtpunkt- og trapesmetoden kan en numerisk likningsløser komme til nytte.)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}