Vecka 7, inför föreläsningen

Denna sista veckan kommer vi att titta på metoder för att lösa ekvationer, beräkna integraler mm numeriskt. Vi såg ju förra veckan att solveset misslyckades ibland att lösa ekvationer symboliskt även då det fanns symboliska lösningar, och för de "flesta" ekvationer så finns det inte ens något sätt att hitta exakta lösningar utan man är hänvisad till att bestämma en approximativ lösning numeriskt. Innan vi tittar på detta så ska vi gå igenom två finesser i Python som kan vara användbara denna veckan (och i framtiden).

Lambda-funktioner

Det finns ett kortfattat och enkelt sätt att definiera funktioner i Python med kommandot Links to an external site.lambda Links to an external site.. Speciellt praktiskt är det när funktionen bara används en enstaka gång som t.ex. argument till en annan funktion vilket vi kommer att göra mycket denna veckan. Om man t.ex. vill definiera andragradsfunktionen f(x)=x²-x-1 så kan man göra så här

f = lambda x : x**2 - x - 1
print(f(3), f(0), f(-3))

vilket ger utskriften

5 -1 11.

Detta är ekvivalent med att göra den vanliga funktionsdefinitionen:

def f(x):
    return x**2 - x - 1

Frivilliga parametrar

Något som ni redan använt mycket i inbyggda funktioner i Python är frivilliga parametrar, dvs parametrar till en funktion som man kan men inte måste ange. Dessa har då ett standardvärde som parametern har om den inte anges. Ett exempel är funktionen plot i modulen matplotlib.pyplot Links to an external site.som har många möjliga parametrar. Dels behöver man inte ange alla de tre grundläggande parametrarna x, y och fmt utan det kan räcka med x och y eller enbart y. Dessutom finns det ett stort antal parametrar som kan anges med sitt namn.

plot(x, y, fmt)
plot(y)                                 # x sätts till range(0, len(y)) och fmt till '-b' (blå linje)
plot(x, y, alpha = 0.1, linewidth = 5)  # Ger bredare genomskinlig blå linje

Man kan när man definierar egna funktioner också ange att vissa parametrar ska vara frivilliga. Dessa anges då efter alla obligatoriska paramterar tillsammans med sitt standardvärde. Här är ett enkelt exempel på en funktion minrot som har en frivillig parameter grad som har 2 som standardvärde.

def minrot(x, grad = 2):
    return x ** (1/grad)

Alla följande anrop av denna funktionen kommer att beräkna kvadratroten ur 5

minrot(5), minrot(5, 2), minrot(5, grad = 2), minrot(x = 5), minrot(x = 5, grad = 2), minrot(grad = 2, x = 5).

Observera att vi här kan välja att anropa funktionen med variablernas namn istället för i en viss ordning och att ordningen då inte spelar någon roll. (Det är långt ifrån så för alla inbyggda funktioner i Python att alla parametrar kan anropas med sitt namn.) Följande anrop kommer istället att beräkna kubikroten ur 5

minrot(5, 3), minrot(5, grad = 3), minrot(x = 5, grad = 3), minrot(grad = 3, x = 5).

Lösa ekvationer numeriskt

Jag kommer i detta avsnittet att referera flitigt till boken som används i matematikkursen på Globala system. Lärarstudenter får antingen konsultera den litteratur man använt tidigare alternativrt ta en titt på wikipedia-sidan Links to an external site. om algoritmer för att hitta rötter till skalära ekvationer (rekommenderas också för er som går Globala system). De viktigaste formlerna finns också i texten nedan.

I kapitel 6 i matteboken finns en (utmärkt) genomgång av grunderna för att lösa ekvationer numeriskt. I avsnitt 6.1 beskrivs hur alla ekvationer kan skrivas på normalform f(x)=0 och på fixpunktsform x=g(x) och hur man kan gå från den ena formen till den andra. I korthet så kan man t.ex. gå från fixpunktsform x=g(x) till normalform genom att subtrahera g(x) från båda leden, f(x)=x-g(x)=0, och från normalform till fixpunktsform genom att multiplicera båda leden med nollskild a(x) och addera x till båda leden

$LaTeX: f(x)=0 \Longleftrightarrow a(x)f(x)=0 \Longleftrightarrow x= x+a(x )f(x) \text{ om } a(x)\neq0$ $f(x)=0 \Longleftrightarrow a(x)f(x)=0 \Longleftrightarrow x= x+a(x )f(x) \text{ om } a(x)\neq0$ .

Bisektionsalgoritmen

I avsnitt 6.2 i boken beskrivs bisektionsalgoritmen (kallas också för intervallhalveringsmetoden) som utgår från normalformen av en ekvation och som är en enkel och stabil metod att lösa ekvationer som garanterat fungerar för alla kontinuerliga funktioner. Idén bygger på att om f(a)f(b)<0, dvs f(a) och f(b) har olika tecken så måste f(x)=0 för något x i intervallet (a,b) om f är kontinuerlig. (Detta är ett specialfall av satsen om mellanliggande värden.) Metoden går ut på att successivt halvera intervallet som innehåller lösningen till ekvationen tills det är så litet som man önskar och då vet man att lösningen ligger i ett mycket litet intervall och att felet är högst halva längden av intervallet om man tar mittpunkten som approximativ lösning.

Första steget är att beräkna f((a+b)/2) och då gäller att antingen är f((a+b)/2)f(a)<0, f((a+b)/2)f(b)<0 eller f((a+b)/2)=0 (Bingo!). Om det första inträffar så fortsätter man med intervallet (a, (a+b)/2). Om det andra inträffar så fortsätter man med intervallet ((a+b)/2, b). Om det tredje inträffar har vi hittat en lösning. Algoritmen finns beskriven med så kallad pseudokod på sidan 198 i matteboken. En av uppgifterna denna vecka är att realisera denna med Python-kod.

En uppenbar nackdel med bisektionsalgoritmen är förstås att man först måste identifiera ett intervall som innehåller en lösning, men det kan man oftast åstadkomma genom att beräkna f(x) i väl valda punkter eller att rita dess graf på lämpliga intervall. En annan begränsning är att den kommer att returnera en enda lösning oavsett hur många lösningar det finns i intervallet (a,b).

Fixpunktsalgoritmen

Avsnitt 6.3 i boken handlar om fixpunktsalgoritmen. Denna utgår från ekvationen på en fixpunktsform x=g(x) och bestämmer en lösning genom att beräkna x_k för följden

$LaTeX: \left\{ \begin{eqnarray*} x_1 &=& \text{gissning} \\ x_k &=&g(x_{k-1}), \text{ för } k>1 \end{eqnarray*} \right.$ $\left\{ \begin{eqnarray*} x_1 &=& \text{gissning} \\ x_k &=&g(x_{k-1}), \text{ för } k>1 \end{eqnarray*} \right.$

Ni har redan gjort ett exempel med denna algoritm i laborationen i vecka 3 då ni beräknade roten ur ett tal n. Roten ur n är en lösning till x²-n=0 och den omskrivning ni använde under vecka 3 var följande

$LaTeX: 0=x^2-n \Longleftrightarrow 0=-\frac1{2x}(x^2-n) \Longleftrightarrow x = x -\frac1{2x}(x^2-n) =x-\frac x2+\frac n{2x}=\frac x2+\frac n{2x} = \frac{x+n/x}{2}=g(x)$ $0=x^2-n \Longleftrightarrow 0=-\frac1{2x}(x^2-n) \Longleftrightarrow x = x -\frac1{2x}(x^2-n) =x-\frac x2+\frac n{2x}=\frac x2+\frac n{2x} = \frac{x+n/x}{2}=g(x)$ .

Om man jämför med den allmänna omskrivningen från normalform till fixpunktsform ovan så är alltså LaTeX: a(x)=-1/2x $a(x)=-1/2x$ . Här är

$LaTeX: g'(x)=\frac12-\frac{n}{2x^2} \text{ så att } g'(\sqrt n)=\frac12-\frac{n}{2\sqrt n^2}=0$ $g'(x)=\frac12-\frac{n}{2x^2} \text{ så att } g'(\sqrt n)=\frac12-\frac{n}{2\sqrt n^2}=0$

vilket betyder att om man är nära så kommer den att konvergera mycket snabbt vilket ni såg när ni körde den i laborationen under vecka 3.

En (kanske lite väl) enkel implementering av fixpunktsalgoritmen i Python kan se ut så här

def fixpunkt(g, x0, tol = 1e-10):
    y = x0          # Startvärde
    x = y + 2 * tol # Ser till så att |x-y| > tol så att den körs minst 1 gång
    while abs(x-y) > tol:
        x = y     # Sätt x till senast beräknade värde
        y = g(x)
    return y

Här är g funktionen i x=g(x), x0 startvärdet för algoritmen och tol en frivillig parameter som anger hur stort avstånd det får högst vara mellan x_n och x_n-1 för att avbryta. Observera speciellt att en av parametrarna, g, vi skickar in i funktionen fixpunkt är själv en funktion. Denna implementering har en stor brist och det är att den riskerar att hålla på för evigt, eftersom det inte är säkert att villkoret i while-loopen någonsin kommer att bli uppfyllt. Man bör ha ett extra villkor där som ser till att loopen avbryts efter ett visst antal iterationer så att man ger upp istället för att fortsätta för evigt.

Ett anrop för att beräkna roten ur n med metoden som beskrevs ovan kan se ut så här om man använder en lambda-funktion

n = float(input('Mata in ett tal: '))
roten = fixpunkt(lambda x : (x + n/x)/2, 1)
print('Roten ur', n, 'är', roten)

eller så här med en "traditionell" funktion

n = float(input('Mata in ett tal: '))
def gn(x):
    return (x + n/x)/2
roten = fixpunkt(gn, 1)
print('Roten ur', n, 'är', roten)

Detta är ett bra exempel på när det är enkelt och elegant att använda en lambda-funktion.

Newtons metod

Ett specialfall av fixpunktsalgoritmen är Newtons metod som beskrivs i avsnitt 6.4 i matteboken. Denna bygger på att man från normalformen f(x)=0 skriver om på fixpunktsform med a(x)=-1/f'(x) så att man får g(x)=x-f(x)/f'(x). Denna kommer att fungera bra om man inte råkar ut för att f'(x) blir litet (f'(x)=0 blir ju katastrof) för då finns risk att den inte kommer att konvergera. Pseudokod för Newtons metod finns på sidan 208 i boken. En av uppgifterna denna veckan kommer vara att implementera en liten variant på denna i Python.

I exemplet ovan med beräkning av roten ur n så hade vi f(x)=x²-n så att f'(x)=2x och därmed är -1/f'(x)=-1/2x så att den omskrivning som användes där var helt enkelt Newtons metod för f(x)=x²-n.

Beräkna integraler

I avsnitt 2.5 i den andra matteboken ("den blå boken") så beskrivs tre grundläggande metoder för att beräkna värdet av bestämda integraler på ett intervall numeriskt. Alternativt kan man titta på wikipedia-sidan om numerisk beräkning av integraler Links to an external site.. Alla metoderna bygger på att man delar upp intervallet i delar och på var och en av dessa delar approximerar man funktionen med en enkel funktion. För mittpunktsmetoden med en konstant funktion, med trapetsmetoden med en linjär funktion och med Simpsons formel med en kvadratisk funktion.

Om vi antar att man delar upp intervallet [a,b] i N lika stora intervall så kommer längden av varje delintervall att vara d=(b-a)/N. Punkterna som ligger i kanterna av de olika delintervallen kommer att vara x₀=a, x₁=a+d, x₂=a+2d, ..., x_N-1=a+(N-1)d=b-d, x_N=a+Nd=b.

Indelning i intervall.png

Med dessa beteckningar så kommer delintervall i att gå från x_i-1till x_i.Formeln för mittpunktsmetoden, som bygger på att man tar funktionen värde mitt i intervallet och multiplicerar med längden av intervallet, blir då

$LaTeX: \int_a^b f(x) dx \approx \sum_{i=1}^N f\left(\frac{x_i+x_{i-1}}{2} \right)\cdot d = d\cdot \sum_{i=1}^N f\left(a+d\left(i-\frac12\right)\right)$ $\int_a^b f(x) dx \approx \sum_{i=1}^N f\left(\frac{x_i+x_{i-1}}{2} \right)\cdot d = d\cdot \sum_{i=1}^N f\left(a+d\left(i-\frac12\right)\right)$ .

För trapetsmetoden tar man istället trapetsen som man får genom att dra en rät linje mellan funktionens värde i start- och slutpunkt. Arean för den trapetsen blir medelvärdet av funktionsvärdena multiplicerat med längden av intervallet dvs.

$LaTeX: \int_a^b f(x) dx \approx \sum_{i=1}^N \frac{f(x_i)+f(x_{i-1})}{2} \cdot d = d\cdot \sum_{i=1}^N \frac{f(a+id)+f(a+(i-1)d)}2 = d\left( \frac{f(a)+f(b)}2 + \sum_{i=1}^{N-1} f(a+id)\right)$ $\int_a^b f(x) dx \approx \sum_{i=1}^N \frac{f(x_i)+f(x_{i-1})}{2} \cdot d = d\cdot \sum_{i=1}^N \frac{f(a+id)+f(a+(i-1)d)}2 = d\left( \frac{f(a)+f(b)}2 + \sum_{i=1}^{N-1} f(a+id)\right)$ .

I den sista omskrivningen utnyttjas att f(a+kd)/2 för k=1,...,N-1 kommer att vara med i summan två gånger då i=k och i=k+1, men bara en gång då k=0 samt k=N.

SciPy

OBS: Detta kommer att gås igenom noggrant på föreläsningen.

Det finns en stor (och växande) modul SciPy Links to an external site. som utnyttjar bl.a. NumPy och matplotlib och som innehåller funktionerna som beskrivits här och mycket, mycket mer. Funktioner för att lösa ekvationer finns i delmodulen scipy.optimize Links to an external site. och funktioner för att beräkna integraler i delmodulen scipy.integrate Links to an external site..

Lösa ekvationer

Den grundläggande funktionen i optimize för att hitta rötter till en (skalär) funktion f, dvs lösa f(x)=0, är root_scalar Links to an external site.. För denna funktion kan man välja vilken metod som ska användas med bisektionsalgoritmen och Newtons metod som två av alternativen. Vi ska se hur man anropar denna för att lösa ekvationen f(x) = 3cos(x)+ln(x)=0. Vi börjar med att plotta funktionen för att få en uppfattning om var det finns rötter. Observera att f bara är definierad för x>0 och när ln(x)>3, dvs x>e³ som är drygt 20, så kan det inte finnas någon rot eftersom |cos(x)|<1. Det räcker alltså att plotta funktionen på intervallet (0,25):

Vi ser att det ser ut att finnas 7 lösningar till ekvationen. Den andra roten ser ut att ligga mellan 1 och 3 så vi testar om bisektionsalgoritmen kan hitta denna och om vi använder root_scalar så ser koden ut så här:

import numpy as np
from scipy import optimize

f = lambda x : 3 * np.cos(x) + np.log(x)
solbisek = optimize.root_scalar(f, bracket = [1, 3], method = 'bisect')   # Sök lösning i intervallet [1,3]
print(solbisek)            # Skriver ut de fem saker som returneras: converged, flag, function_calls, iterations, root
if solbisek.converged:     # Om metoden konvergerar så skriv ut roten och räkna ut funktionens värde i denna
    print("\nRot: ", solbisek.root, "Värde för funktionen i punkten:",  f(solbisek.root))
else:                      # Om metoden inte konvergerar så skriv ut felmeddelandet
    print(solbisek.flag)

vilket ger utskriften

      converged: True
           flag: 'converged'
 function_calls: 42
     iterations: 40
           root: 1.7604555043781147

Rot:  1.7604555043781147 Värde för funktionen i punkten: -1.6249224188413791e-12

Observera att root_scalar förutom en rot också returnerar om den konvergerade (vilket är mycket viktig info) , hur många iterationer som behövdes samt hur många gånger som f anropats. Om man väljer ett intervall där funktionen inte har olika värden i ändpunkterna, t.ex. sätter bracket = [2,3] så får man felmeddelandet "ValueError: f(a) and f(b) must have different signs"

Om man istället vill använda Newtons metod så behöver man också derivatan av funktionen. Koden ser ut så här om man tar x₀=1 som startvärde:

import numpy as np
from scipy import optimize

f = lambda x : 3 * np.cos(x) + np.log(x)
fp = lambda x : - 3 * np.sin(x) + 1/x
solnewton = optimize.root_scalar(f, x0 = 1, fprime = fp, method = 'newton')
print(solnewton)
if solnewton.converged:
 print("\nRot: ", solnewton.root, "Värde för funktionen i punkten:", f(solnewton.root))
else:
 print(solnewton.flag)

vilket ger utskriften

      converged: True
           flag: 'converged'
 function_calls: 10
     iterations: 5
           root: 1.7604555043774315

Rot:  1.7604555043774315 Värde för funktionen i punkten: -1.1102230246251565e-16

Observera att Newtons metod behövde mycket färre iterationer och hittade en bättre approximation (funktionsvärde mycket mindre). När den fungerar är Newtons metod mycket snabbare än bisektionsalgoritmen. Observera dock att Newtons metod är känslig för val av startpunkt. Om vi stället väljer x₀=0.5 och kör exakt samma kod så kommer den att misslyckas med att konvergera och ett felmeddelande avslöjar att det blir ett negativt x under kalkylen och utskriften blir nu

      converged: False
           flag: 'convergence error'
 function_calls: 100
     iterations: 50
           root: nan

convergence error
C:\XXX.py:79: RuntimeWarning: invalid value encountered in log
  f = lambda x : 3 * np.cos(x) + np.log(x)

Beräkna integraler

Den generella metoden i scipy.integrate för att beräkna bestämda integraler (med en variabel) är quad. Denna använder olika metoder beroende på funktion och intervall. Anropet är enkelt och för att beräkna

$LaTeX: \int_{0}^{\pi} \frac{\sin(x)}{x} dx$ $\int_{0}^{\pi} \frac{\sin(x)}{x} dx$

är anropet

import numpy as np
from scipy import integrate

int1 = integrate.quad(lambda x : np.sin(x)/x, 0, np.pi)
print(int1)

som ger utskriften

(1.851937051982466, 2.0560631552673694e-14).

Här är det första talet det beräknade approximativa värdet på integralen och det andra en övre gräns på felet.

Man kan också beräkna generaliserade integraler med gränser som är plus och eller minus oändligheten så att t.ex. beräkna

$LaTeX: \int_{-\infty}^{\infty} e^{-x^2} dx$ $\int_{-\infty}^{\infty} e^{-x^2} dx$

genom följande kod

import numpy as np
from scipy import integrate

int1 = integrate.quad(lambda x : np.exp(-x**2), -np.inf, np.inf)
print(int1)

vilket ger utskriften

(1.7724538509055159, 1.4202636781830878e-08).

(Om man kvadrerar det approximativa värdet på integralen så får man $LaTeX: \pi$ $\pi$ med 14 korrekta decimaler. Funktionen $LaTeX: e^{-x^2}$ $e^{-x^2}$ saknar elementär primitiv funktion, men värdet på just denna bestämda integral går trots det att bestämma exakt till $LaTeX: \sqrt{\pi}$ $\sqrt{\pi}$ .)