1.1. Ber¨ akning av integraler

(1)

Kapitel 1. Numeriska metoder

Detta är andra delen av kursen i vetenskapliga beräkningar, där vi till en början kommer att bekanta oss med endel numeriska metoder, som inte ingick i den första delen. Beräkningarna kommer som förut att utföras med Matlab.

(2)

1.1. Ber¨ akning av integraler

I fysiken behöver man ofta beräkna integraler av funktioner numeriskt. För att beräkna integraler av elementära funktioner finns det standardmetoder som är kända fr˚an analysen, men i fysiken stöter man ofta p˚a funktioner, som inte kan integreras med välkända metoder.

Ett exempel h¨arp˚a ¨ar Gauss’ felfunktion

erf(x) = 2

√π Z x

0

e^−t2dt

som är nära besläktad med den kumulativa normalfördelningen som används i statistiken. Denna integral kan dock inte beräknas med elementära metoder eftersom den primitiva funktionen till e^−t2 inte kan uttryckas med elementära funktioner. Naturligtvis kunde man kalla funktionen för erf, men hjälper oss inte att beräkna värden av den. För den skull behöver man numeriska metoder. Redan i antiken fanns det praktiskt behov av att beräkna ytor och volymer. Arkimedes uppfann en metod för s˚adana beräkningar som i hög grad p˚aminner om beräkning av integraler.

Den typ av integraler, som man normalt st¨oter p˚a i fysiken, kallas i matematiken f¨or en Riemann-integral.

(3)

Riemann-integralen baserar sig p˚a Jordanm˚attet, och definieras som gr¨ansv¨ardet av en Riemannsumma, som vi skall se nedan.

Om vi för enkelhetens skull antar att vi önskar beräkna integralen över funktionen f (x) mellan gränserna x = a och x = b, s˚a kan vi tänka oss detta intervall för enkelhetens skull delat i n lika stora delintervall av längden h = (b − a)/n. Delningspunkterna kan d˚a betecknas x_k = a + kh, (k = 0, 1, . . . , n), dvs h = x_k − x_k−1. En Riemannsumma kan d˚a definieras som

S =

n

X

k=1

f (c_k)(x_k − x_k−1), c_k ∈ [x_k−1, x_k]

Varje term i denna summa är allts˚a ytan av en rektangel med höjden f (c_k) och bredden h, och summan kan uppfattas som en approximation till en yta som begränsas av x-axeln, funktionens ordinator i intervallets ändpunkter, och funktionskurvan mellan ordinatorna. En funktion sägs vara Riemannintegrerbar, ifall Riemannsumman konvergerar mot ett bestämt gränsvärde, d˚a delningen tätnar.

Ofta brukar man välja ck s˚a, att det antingen är funktionens maximivärde (supremum) Mk inom intervallet [x_k, x_k−1] eller minimivärde (infimum) mk inom detta intervall. Vi f˚ar d˚a en övre Darbouxsumma U = Pn

k=1M_k(x_k−xk−1), respektive en undre Darbouxsumma L = Pn

k=1m_k(x_k−xk−1). D˚a delningen tätnar, kommer den övre summan att uppifr˚an närma sig ett gränsvärde, som kallas Darboux’ övre integral, och den undre summan närmar sig nedifr˚an ett gränsvärde, som kallas Darboux’ undre integral.

(4)

Om b˚ada integralerna sammanfaller, s˚a kallas det gemensamma gr¨ansv¨ardet Darboux’ integral, som d˚a sammanfaller med Riemanns integral.

Riemannsummor kan ocks˚a använda för att numeriskt beräkna en integral. Som ett exempel skall vi ta R 1

0 e^xdx. F¨ordenskull delar vi upp intervallet [0, 1] p˚a n delintervall av l¨angden h = 1/n. Den undre Darbouxsumman blir d˚a

L =

n

X

k=1

m_k(x_k − x_k−1) =

n

X

k=1

e^(k−1)/n1

n = 1 n

n

X

k=1

(e^1/n)^k−1

=1 n

(e^1/n)ⁿ − 1

e^1/n − 1 = (e − 1) 1/n

e^1/n − 1, (1)

eftersom funktionen är monotont växande och serien är geometrisk. Om delningen förtätas gränslöst, dvs man l˚ater n → ∞, s˚a blir Darbouxsummans gränsvärde e − 1, vilket är integralens värde.

(5)

1.2. Kvadraturregler

I allmänhet är det lättare att beräkna integraler numeriskt med hjälp av s.k. kvadraturregler än med Riemannsummor. I detta fall använder man en vägd summa av funktionsvärden (som kallas noder) inom det givna intervallet för att approximera integralen. Vikterna motsvarar d˚a intervallängderna i Riemannsumman.

En m punkters kvadraturregel Q f¨or den best¨amda integralen I =

Z b a

f (x)dx

¨ar en approximation, som har formen

Q = (b − a)

m

X

k=1

w_kf (x_k),

där xk kallas för abskissor och w_k är vikterna. Regeln är definierad, d˚a abskissor och vikter är angivna.

Effektiviteten för en kvadraturregel mäts med antalet funktionsberäkningar som krävs. En sex-punkters formel är s˚aledes dubbelt dyrare än en tre-punkters formel.

(6)

Vi skall börja med att studera formler av Newton–Cotes typ. De kännetecknas av att man integrerar en polynomisk interpolant p(x) för integranden f (x). Detta betyder, att av p(x) ≈ f (x) följer

Z b a

f (x)dx ≈ Z b

a

p(x)dx.

Newton–Cotes kvadraturregler f˚ar man genom att integrera polynom, som interpolerar integranden i punkter p˚a lika avst˚and fr˚an varandra. Newton–Cotes m-punktsformel definieras av

Q^(m)_{N C} = Z b

a

p_m−1(x)dx,

d¨ar p_m−1(x) interpolerar f (x) i punkterna

x_i = a + i − 1

m − 1(b − a), i = 1 : m.

(7)

Om m = 2, f˚ar vi trapetsregeln Q⁽²⁾_{N C} =

Z b a

f (a) + f (b) − f (a)

b − a (x − a)

dx

= (b − a)₁

2f (a) + ¹₂f (b) . (2)

Om m = 3 och c = (a + b)/2, s˚a f˚ar vi Simpsons regel:

Q⁽³⁾_{N C} = Z b

a

"

f (a) + f (c) − f (a)

c − a (x − a) +

f (b)−f (c)

b−c − f (c)−f (a) c−a

b − a (x − a)(x − c)

# dx

= b − a 6

f (a) + 4f a + b 2

+ f (b)

. (3)

Av dessa exempel framg˚ar, att kvadraturformeln kan uttryckas som en lineärkombination av funktionsvärden p˚a jämnt avst˚and fr˚an varandra innanför integrationsintervallet.

(8)

Ett enkelt sätt att finna vikterna i en dylik kvadraturformel f˚ar man genom att fordra, att den skall stämma exakt för x–potenserna xⁱ, i = 0 : m − 1. Om vi t.ex. vill beräkna koefficienterna i fyrpunktsformeln:

Z 1 0

f (x)dx ≈ af (0) + bf 1 3

+ cf 2 3

+ df (1),

s˚a st¨aller vi upp ekvationssystemet

Z 1 0

dx = 1 = a + b + c + d Z 1

0

xdx = 1

2 = 0 + b

3 + 2c

3 + d Z 1

0

x²dx = 1

3 = 0 + b

9 + 4c

9 + d Z 1

0

x³dx = 1

4 = 0 + b

27 + 8c

27 + d. (4)

(9)

Genom att l¨osa detta ekvationssystem finner man formeln Z 1

0

f (x)dx ≈ 1 8

f (0) + 3f 1 3

+ 3f 2 3

+ f (1)

.

F¨or enkelhetens skull kan vi skriva ett litet MATLAB–program, som tabulerar de vanligaste Newton–Cotes–

koefficienterna:

function w=wnc(m);

%

% Invariabel:

% m heltal mellan 2 och 8

%

% Utvariabel:

% w kolumnvektor med m element som inneh˚aller vikterna

% f¨or m-punktsregeln.

%

if m==2

w = [1 1]’/2;

elseif m==3

w = [1 4 1]’/6;

elseif m==4

w = [1 3 3 1]’/8;

(10)

elseif m==5

w = [7 32 12 32 7]’/90;

elseif m==6

w=[19 75 50 50 75 19]’/288;

elseif m==7

w=[41 216 27 272 27 216 41]’/840;

elseif m==8

w=[751 3577 1323 2989 2989 1323 3577 751]’/17280;

end

Som vi ser, är viktsvektorerna symmetriska kring intervallets mittpunkt, s˚a att w(1 : m) = w(m : −1 : 1). Den numeriska integrationen kan sedan utföras med hjälp av formeln

Q^(m)_{N C} = (b − a)

m

X

i=1

w_if_i = (b − a)[w₁ · · · w_m]





f (x₁) ...

f (x_m)





som visar att kvadraturformeln kan beräknas som skalärprodukten av viktsvektorn och en vektor inneh˚al- lande funktionsvärden. En MATLAB-rutin, som beräknar en integral med hjälp av Newton–Cotes kvadratur, kan allts˚a skrivas

(11)

function I = qnc(fname,a,b,m)

%

% Invariabler:

% fname : str¨ang som inneh˚allet namnet p˚a en funktion f(x) som ¨ar

% definierad ¨over intervallet [a,b]. f ger till resultat

% en kolumnvektor, om x ¨ar en kolumnvektor.

% a,b : reella skal¨arer.

% m : ett heltal mellan 2 och 8.

% Utvariabel:

% I : en m-punkts Newton-Cotes approximation f¨or integralen

% av f(x) mellan a och b.

%

w = wnc(m);

x = linspace(a,b,m)’

f = feval(fname,x);

I = (b-a)*(w’*f);

Följande program kan användas för att testa funktionen:

(12)

% Skriptfil: testqnc

%

% Testar Newton-Cotes-reglerna.

while input (’Annat exempel? (1=ja, 0=nej). ’)

fname = input(’Ange namnet p˚a integranden mellan apostrofer:’);

a = input(’Ange v¨anstra ¨andpunkten: ’);

b = input(’Ange h¨ogra ¨andpunkten: ’);

s = [’QNC(’ fname sprintf(’,%6.3f,%6.3f,m)’,a,b)];

clc

disp([’ m ’ s])

disp(’ ’) for m=2:8

I = qnc(fname,a,b,m);

disp(sprintf(’ %2.0f %20.16f’,m,I)) end

end

Om programmet till¨ampas p˚a sinusfunktionen inom intervallet [0, π/2] f˚as

(13)

m QNC(sin, 0.000, 1.571,m )

2 0.7853981633974483

3 1.0022798774922104

4 1.0010049233142790

5 0.9999915654729928

6 0.9999952613861668

7 1.0000000258372352

8 1.0000000158229038

Som vi ser, ¨okar noggrannheten med antalet punkter i formeln. Man kan visa, att felet i Simpsons regel (m = 3) best¨ams ur formeln

Z b a

f (x)dx − Q⁽³⁾_{N C}

≤ (b − a)⁵

2880 M₄,

där M4 betecknar övre gränsen av fjärde derivatan |f⁽⁴⁾(x)| inom intervallet [a, b]. I det allmänna fallet gäller formeln

Z b a

f (x)dx = Q^(m)_{N C} + c_mf^(d+1)(θ) b − a m − 1

^d+2 ,

(14)

d¨ar cm ¨ar en liten konstant, θ ∈ [a, b] och

d =

(m − 1 , om m är jämnt m , om m är udda.

Om vi känner till en övre gräns för derivatan, t.ex. att |f^(d+1)(x)| ≤ M_d+1 inom intervallet [a, b], s˚a

¨ar

Q^(m)_{N C} − Z b

a

f (x)dx

≤ |c_m|M_d+1 b − a m − 1

^d+2 .

Här är en MATLAB–funktion som räknar ut kvadraturfelet d˚a känner en uppskattning av felet i derivatan:

function error = ncerr(a,b,m,DerBound)

%

% Invariabler:

% a,b reella tal som satisfierar a<=b

% m heltal som satisfierar 2<=m<=8 och

% DerBound en övre gräns för den d:te derivatan av

% en funktion f(x) definierad p˚a [a,b]. Talet

% d = m om m är udda, och m-1 om m är jämnt.

%

% Utvariabel:

% error en övre gräns för absoluta felet i en Newton-

(15)

% Cotes m-punktsregel, d˚a den till¨ampas p˚a

% integralen av f(x) fr˚an a till b.

%

if m==2

d=1; c = -1/12;

elseif m==3

d=3; c = -1/90;

elseif m==4

d=3; c = -3/80;

elseif m==5

d=5; c = -8/945;

elseif m==6

d=5; c = -275/12096;

elseif m==7

d=7; c = -9/1400;

elseif m==8

d=7; c= -8183/518400;

end;

error = abs( c*DerBound*((b-a)/(m-1))^(d+2) );

Rutinen visar, att en (m − 1)-punktsregel vanligen är lika bra som motsvarande m-punktsregel. Vi ser detta tydligare genom att köra följande testprogram:

(16)

% Skriptfil: shqncerr

%

% Program f¨or att studera felet i en Newton-Cotes regel.

clc

disp(’Enkelt fall: integralen fr˚an 0 till pi/2 av sin(x)’) disp(’ ’)

disp(’Antag att DerBound = 1.’) disp(’ ’)

disp(’ m QNC(m) fel felgr¨ans’)

disp(’ ’) for m=2:8

numI = qnc(’sin’,0,pi/2,m);

err = abs(numI-1);

errBound = errnc(0,pi/2,m,1);

s = sprintf(’%20.16f %10.3e %10.3e’,numI,err,errBound);

disp([ sprintf(’ %2.0f ’,m) s]) end

Resultatet f¨or integralen av sinus mellan 0 och π/2 ser ut s˚a h¨ar:

(17)

m QNC(m) fel felgr¨ans 2 0.7853981633974483 2.146e-001 3.230e-001 3 1.0022798774922104 2.280e-003 3.321e-003 4 1.0010049233142790 1.005e-003 1.476e-003 5 0.9999915654729928 8.435e-006 1.219e-005 6 0.9999952613861668 4.739e-006 6.867e-006 7 1.0000000258372352 2.584e-008 3.714e-008 8 1.0000000158229038 1.582e-008 2.277e-008

Felet i kvadraturegeln är beroende av längden p˚a intervallet, som vi sett. Vi minskar därför p˚a felet, om vi kan minska p˚a intervallängden. Detta ˚astadkoms genom att dela intervallet [a, b] i delintervall. Om vi antar att a = z₁ < z₂ < · · · < z_n+1 = b, s˚a kan integralen beräknas som en summa av integraler, tagna över delintervallen:

Z b a

f (x)dx =

n

X

i=1

Z _zi+1

zi

f (x)dx.

(18)

Om vi till¨ampar Q^(m)_{N C} p˚a var och en av delintegralerna, s˚a f˚ar vi en sammansatt kvadraturregel. Om t.ex.

∆_i = z_i+1 − zi och z

i+1 2

= (z_i + z_i+1)/2, i = 1 : m, s˚a g¨aller

Q =

n

X

i=1

∆_i 6

f (z_i) + 4f (z

i+1 2

) + f (z_i+1)

.

Om z är en MATLAB–vektor, som inneh˚aller delningspunkterna i [a, b], och fname är en sträng, som inneh˚aller namnet p˚a en funktion, s˚a kommer MATLAB–programfragmentet

I = 0

for i=1:length(z)-1

I = I + qnc(fname,z(i),z(i+1),m);

end

att ge I värdet av en sammansatt m–punkters Newton-Cotes approximation till integralen, som baserar sig p˚a delningspunkterna i z. Man kan ocks˚a välja delningspunkterna p˚a ett mera automatiskt sätt. Detta utnyttjas vid adaptiv kvadratur, där delningspunkterna väljs tätare, d˚a funktionen varierar snabbare, och glesare, d˚a funktionen varierar l˚angsammare (se t.ex. MATLAB–procedurerna quad och quadl).

(19)

I Gauss kvadratur väljs abskissorna p˚a ett s˚adant sätt, att regeln blir exakt för polynom av största möjliga gradtal. Som ett exempel skall vi bestämma vikterna och abskissorna s˚a att formeln

w₁f (x₁) + w₂f (x₂) = Z 1

−1

f (x)dx

är exakt för polynom av gradtalet 3 eller mindre. Genom att kräva att regeln skall stämma exakt för de fyra funktionerna 1, x, x², x³ f˚ar vi fyra ekvationer med fyra obekanta (jfr v˚ar beräkning av vikterna i en trepunkters Newton-Cotesregel):

w₁ + w₂ = 2 w₁x₁ + w₂x₂ = 0 w₁x²₁ + w₂x²₂ = 2 3

w₁x³₁ + w₂x³₂ = 0. (5)

Detta ekvationssystem kan lösas genom att man multiplicerar den andra ekvationen medx²₁ och subtraherar den fr˚an den fjärde ekvationen. Vi f˚ar d˚a w₂x₂(x₁² − x²₂) = 0, varav följer x₂ = −x₁. Av den andra ekvationen följer d˚a w₁ = w₂, vilket i kombination med den första ekvationen leder till w₁ = w₂ = 1.

(20)

Den tredje ekvationen ger x²₁ = 1/3, varav f¨oljer x₁ = −1/√

3 och x₂ = 1/√

3. Vi finner slutligen Z 1

−1

f (x)dx ≈ f (−1/√

3) + f (1/√ 3), vilket kallas tv˚apunkters Gauss–Legendre regeln.

En m-punkters Gauss–Legendre regel har formen

Q^(m)_GL = w₁f (x₁) + · · · + w_mf (x_m),

där vikterna och abskissorna valts s˚a, att regeln är exakt för polynom av gradtalet 2m − 1.

H¨arnedan visas en MATLAB–rutin, som ger vikter och abskissor f¨or upptill 4 punkters Gauss–Legendre regler:

function [w,x] = wgl(m);

%

% Invariabel:

% m heltal som uppfyller villkoret 2 <= m <= 4

%

% Utvariabler:

% w kolumnvektor med m element som inneh˚aller vikterna f¨or

(21)

% en m-punkters Gauss-Legendre regel.

% x en kolumnvektor med m element som inneh˚aller abskissorna

% f¨or en m-punkters Gauss-Legendre regel.

%

w = ones(m,1);

x = ones(m,1);

if m==2

w(1) = 1.000000000000000; w(2) = w(1);

x(1) = -0.577350269189626; x(2) = -x(1);

elseif m==3

w(1) = 0.555555555555558; w(3) = w(1);

w(2) = 0.888888888888889;

x(1) = -0.774596669241483; x(3) = -x(1);

x(2) = 0.000000000000000;

elseif m==4

w(1) = 0.347854845137454; w(4) = w(1);

w(2) = 0.652145154862546; w(3) = w(2);

x(1) = -0.861136311594053; x(4) = -x(1);

x(2) = -0.339981043584856; x(3) = -x(2);

end;