Lösning av ekvationer f(x) = 0
Max och Min
för f(x) på ett intervall [a, b]
Numeriska metoder har fått en ökad betydelse för ingenjörsämnena under senare år. Anledningen är givetvis den snabba utvecklingen inom datorområdet, från persondatorer till superdatorer.
Vårt mål är att utgående från matematiska modeller (funktioner och ekvationer) ta fram algoritmer som vi kan programmera (i C#). Algoritmerna ger en lösning på problemet i form av tal.
Praktisk ingenjörsmatematik slutar förr eller senare med numeriska resultat. Ingenjörsstudenter måste därför lära sig ett antal standardmetoder för att lösa numeriska problem. Det är till och med så att många praktiska problem endast kan lösas numeriskt. Typiska exempel är lösningen av många ekvationer och beräkningen av de flesta integraler.
Det normala är att det inte går att uttrycka rötterna till en ekvation f(x) = 0 i sluten form. Approximativa metoder får därför tillgripas. Dessa utnyttjar i allmänhet iteration.
Det är lämpligt att inleda lösandet av ekvationsproblem med att skaffa sig en grov skattning till den sökta roten. Flertalet förfinade numeriska metoder för ekvationslösning måste utgå från en hyfsad skattning för att konvergera mot roten.
Vi skall för att åskådliggöra de olika metoderna använda följande problem som exempel:
Lös
ekvationen:
Grafisk teknik
Att helt enkelt rita upp f(x) och avläsa
nollställena faller sig ganska naturligt. Även ett enkelt problem kan
emellertid förenklas ytterligare med en mindre omformulering: Att rita upp
grafen till x - e-x är inte
särskilt lätt. Betydligt lättare är det att skriva ekvationen på formen
x = e-x och sedan
i samma figur rita graferna till höger- resp vänsterled. Såväl hur den räta
linjen y = x och hur exponentialkurvan löper känner vi väl till och kan skissa
utan några stora hjälpmedel. Där de båda kurvorna skär varandra avläser man
rötterna till den givna ekvationen. Med den här tekniken erhåller vi alltså
direkt skattningen 0,5 á 0,6 till den sökta roten. Frågor som huruvida det
finns fler rötter till ekvationen kan i allmänhet lätt besvaras sedan man ritat
lite.
Trial and error
Om man önskar skaffa sig en något bättre skattning än den man erhåller grafiskt kan man upprätta en tabell. Med denna arbetar man tills man har nått en noggrannhet som man nöjer sig med i denna inledande fas.
Med dagens effektiva räknehjälpmedel kan det tyckas naturligt att göra något liknande en sådan tabell mha datorn. Man kan använda sig av en trial-and-error-teknik, dvs beräkna f(x) för några x-värden, och hålla på tills man har lokaliserat roten med önskvärd noggrannhet. Att man av en slump skulle finna roten till ekvationen med denna teknik tycks uteslutet, åtminstone torde inte detta ske inom överskådlig tid. Det är inte heller tanken med denna metod. Tanken är att man skall dra slutsatsen att det finns åtminstone en rot mellan x1 och x2 .
Villkoret är att
Det är då naturligt att undersöka hur det förhåller sig mitt emellan dessa punkter. Man fortsätter naturligt på den inslagna vägen, dvs undersöker tecknet på f(x) i mittpunkten på det intervall i vilket man spärrat in roten: Utvecklad bör den här tekniken leda till att roten trängs in i allt trängre intervall. Att successivt just halvera dessa intervall är det naturliga när man endast beaktar tecknet på f(x) i intervallets ändpunkter.
Intervallhalveringsmetoden
Trial and error-tekniken har därmed förvandlats till en systematisk intervallhalverings-metod. Intervallhalveringsmetoden innebär just att successivt halvera det intervall i vilket en rot kan ligga. Till information, om vilket man ska välja av de två intervall man får vid en halvering, använder man endast tecknet på funktionen vars nollställe man söker. Startar man med ett intervall [a, b] sådant att f(a)·f(b) < 0 konvergerar metoden nödvändigtvis mot ett nollställe till f(x). Eftersom man vinner endast en faktor 1/2 i noggrannhet för varje ny undersökning är metoden lite väl långsam för att användas till att bestämma nollställen noggrant. Startar vi t ex med ett intervall av längden 0,5 och önskar roten till f(x) = 0 bestämd med 6 korrekta decimaler, alltså med ett absolutfel som är högst 0,5·10-6, erfordras 20 halveringar innan intervallets längd underskrider 0,5·10-6.
Sekantmetoden
Det är rimligt att tro att om man förutom till tecknet också tar hänsyn till värdet av f(x) i ändpunkterna på det intervall man för tillfället betraktar, så skulle man kunna åstadkomma en effektivare metod. I figuren nedan är x0 och x1 två x-värden som man med trial and error funnit ligga något sånär nära roten a till ekvationen f(x) = 0.
Sekanten genom punkterna (x0, f(x0)) och (x1, f(x1)) skär x-axeln i x2. Räta linjens ekvation ger:
Tanken är att x2 ska skatta roten bättre än vad x0 och x1 gör. Vi får alltså en iterativ metod:
Detta är den s k sekantmetoden för ekvationslösning.
Iterationerna avbryts när värdena inte längre ändras, sett med den noggrannhet man räknar med. Sekantmetoden är avsevärt snabbare än intervallhalveringsmetoden.
Newton-Raphsons metod
Man kan i sekantmetoden formellt låta xn-1 gå mot xn och får därvid
Vi ser att xn-1 är skärningen
mellan x-axeln och tangenten i (x, f(xn)) till
y = f(x).
Den iterationsmetod som nu skisserats är Newton-Raphsons metod.
Metoden har en stor nackdel: Vi måste känna till derivatan!
Sammanfattning
Gör så här för att finna en rot till f(x) = 0:
Ø Använd Trial and error för att finna a och b, sådana att f(a)·f(b) < 0.
Ø Använd sedan Sekantmetoden för att finna roten med önskad noggrannhet.
Som bekant är den analytiska lösningen:
Analytisk lösning
Ø Derivera f(x)
Ø Lös ekvationen 
Ø Rötterna till ekvationen ger Max och Min.
Ø Om det är Max, eller Min bestäms av tecknet hos andraderivatan i rötterna till ekvationen.
Ø Kontrollera f(a) och f(b)
Numerisk lösning
Ø Dela in [a, b] i n (n = "stort") delintervall
Ø Beräkna f(x) i alla punkterna
Ø Det största värdet på f(x) är då Max och vice versa.
En C#-funktion för att beräkna Max kan se ut så här:
static public double MaxOnIntervall(FunctionOneParameter
func, double a, double
b, out double
xmax)
{
const
int n = 1000;
double
max = func(a);
xmax = a;
double
step = (b-a)/n;
double
x = a;
while
(x<= b)
{
x += step;
if
(func(x) > max)
{
xmax = x;
max = func(x);
}
}
return
max;
}
Derivatans definition ger direkt
dvs om vi bildar den sista kvoten och sätter h till ett "litet" tal har vi en approximation till derivatan.
Exempel på en C#-funktion:
static public double Derivative(FunctionOneParameter
func, double x, double
eps)
{
return
(func(x+eps) - func(x-eps))/(2*eps);
}
Kan h vara hur litet som helst? Svar: Nej!
Exempel:
Beräkna derivatan för x3 - x i punkten 4.
Analytisk lösning: 3·42 - 1 = 47
Körning av funktionen ovan med minskande h-värden ger:
h d(x^3-x)/dx
1.0000000000E+00 4.8000000000E+01
1.0000000000E-01 4.7010000000E+01
1.0000000000E-02 4.7000100004E+01
1.0000000000E-03 4.7000000952E+01
1.0000000000E-04 4.6999999904E+01
1.0000000000E-05 4.7000002814E+01
1.0000000000E-06 4.7000008635E+01
1.0000000000E-07 4.7000357881E+01
1.0000000000E-08 4.6993955038E+01
1.0000000000E-09 4.6944478527E+01
1.0000000000E-10 4.6857167035E+01
1.0000000000E-11 4.3655745685E+01
1.0000000000E-12 0.0000000000E+00
Man kanske tror att h bör väljas så litet som möjligt, men det finns flera anledningar att nöja sig med ett h som inte ens tillnärmelsevis är ungefär noll.
En anledning är kancellation av termer i täljaren: När h är mycket litet är de båda funktionsvärdena nästan lika. Om de beräknats med, säg, 6 decimaler varav de 5 första är lika kommer efter subtraktionen endast en decimal att återstå, varför den relativa noggrannheten i täljaren blir mycket dålig.
En annan anledning är att om avrundningsfelen i täljarens funktionsvärden divideras med ett litet tal 2h, så blir deras inflytande på resultatet högst betydande: Om, som vi just antog, funktionsvärdena i täljaren beräknats med 6 korrekta decimaler blir täljarens felgräns 10-6. Om nu h = 10-5 bidrar avrundningsfelen med felet 10-6/(2·10-5) = 0,05 i täljaren, vilket är ett förhållandevis stort fel.
En för nybörjaren i det närmaste oöverskådlig mängd metoder för numerisk integralberäkning finns utvecklade. Man kommer dock förvånansvärt långt med följande enkla metoder.
Approximation med Rektanglar
Värdet av
kan skattas med ytan
av rektanglar enligt figuren nedan:
Om vi låter antalet intervall n gå mot oändligheten kommer den sammanlagda ytan av rektanglarna att närma sig integralvärdet.
Om vi vill "stänga in" integralens värde, bildar vi s k under- och översummor av rektanglar.
Ett bättre sätt att approximera integralen är att approximera den med parallelltrapetser.
Trapetsformeln
Om ytan skattas med parallelltrapetser får vi:
Om denna approximation görs i en följd av intervall, se figuren, får man, om vi betecknar steget med h:
Denna formel för approximation av en bestämd integral kallas trapetsformeln.