f(a + h) = f(a) + f (a)h + f (θ) 2 h2, θ [a, a + h]. = f(a+h) f(a)

(1)

Vi skall nu se, hur man kan beräkna numeriska derivator. Antag att vi vill beräkna derivatan av f (x) i en punkt x = a, och att dess Taylor–utveckling kring denna punkt är

f (a + h) = f (a) + f⁰(a)h + f⁰⁰(θ)

2 h², θ ∈ [a, a + h].

Uttrycket D_h = f (a+h)−f (a)

h kommer därför ge allt noggrannare approximationer för derivatan d˚a h minskar, eftersom D_h = f⁰(a) + f⁰⁰(θ)h/2. Här är ett vektoriserat program med vilket man kan studera felet d˚a man tillämpar formeln p˚a f (x) = sin x för olika h-värden (\n anger ny rad):

a = 1; h = logspace(-1,-16,16);

Dh = (sin(a+h) - sin(a))./h;

fel = abs(Dh - cos(a)); tab = [h;fel];

disp(’ h fel’),fprintf(’ %8.1e %18.15f \n’,tab)

Resultaten lagrades i en matris tab, beh¨andigt utskriven en kolumn per rad med kommandot fprintf :

(2)

h fel

1.0e-01 0.042938553332751 1.0e-02 0.004216324856271 1.0e-03 0.000420825507813 1.0e-04 0.000042074449519 1.0e-05 0.000004207362275 1.0e-06 0.000000420746809 1.0e-07 0.000000041827691 1.0e-08 0.000000002969885 1.0e-09 0.000000052541266 1.0e-10 0.000000058481036 1.0e-11 0.000001168704061 1.0e-12 0.000043240216924 1.0e-13 0.000733915900314 1.0e-14 0.003706976198187 1.0e-15 0.014809206444439 1.0e-16 0.540302305868140

Kommandot fprintf kan ocks˚a användas för att skriva data till en fil (se help fprintf). Vi ser att felet först minskar, men sedan igen börjar öka.

(3)

Varje fel i täljaren multipliceras med en faktor 1/h, s˚a att vi istället för det matematiska felet |D_h − f⁰(a)| ≤ ^h₂|f⁰⁰(θ)| i verkligheten f˚ar felet

|D_h − f⁰(a)| ≈ h

2|f⁰⁰(θ)| + 2eps h , där eps är maskinprecisionen. Högra membrum minimeras, om h = 2p

eps/|f⁰⁰(θ)|.

Vi kan nu konstruera en ny rutin för beräkning av derivator. Om andra derivatan är uppifr˚an begränsad:

|f⁰⁰(x)| ≤ M₂, s˚a gäller följande olikhet för det matematiska felet: |Dh − f⁰(a)| ≤ ^M2₂ h. Om det absoluta felet vid en funktionsberäkning är mindre än δ, s˚a kan det totala felet vid den numeriska beräkningen av funktionsderivatan approximeras med

fel(D(h)) = M₂h

2 + 2δ h .

Detta fel minimeras om h = h_min = 2pδ/M2, vilket ger uttrycket fel(D(h_min)) = M₂pδ/M2+

√ δ

δ/M2 = 2√

δM₂. Vi kan d˚a modifiera programmet p˚a f¨oljande s¨att:

(4)

function d = deriv(fname,a,delta,M2);

%

% Invariabler:

% fname en str¨ang som inneh˚aller funktionens namn.

% a: det värde av x för vilket f’(x) beräknas.

% delta: absoluta felet vid funktionsber¨akningen.

% M2: en uppskattning av andra derivatans storlek n¨ara a.

%

% Utvariabler:

% d: en approximation f¨or f’(a).

% err: en uppskattning av felet i d.

%

hopt = 2*sqrt(delta/M2);

d = (feval(fname,a+hopt) - feval(fname,a))/hopt;

Här har använts funktionen feval, vars första argument är en sträng, som inneh˚aller funktionsnamnet, och de övriga argumenten är funktionens egna argument.

Programmet förutsätter att man kan uppskatta det absoluta felet δ och andra derivatans övre gräns M2, men man kan använda konstanta värden av dessa parametrar för att underlätta programmets användning.

(5)

5.4. Interpolation

D˚a vi vill approximera data, utg˚ar vi fr˚an ett antal punkter (x₁, y₁), (x₂, y₂), . . ., (x_n, y_n) och försöker konstruera en funktion f (x), som beskriver punkterna s˚a bra som möjligt. Typen av funktion bestäms av punkterna. Om de t.ex. oscillerar, s˚a är det lämpligt att använda sig av trigonometriska funktioner för att approximera dem. I m˚anga fall kan ett polynom av l˚agt gradtal vara lämpligt. Minsta kvadratmetoden, som vi redan använt, kan användas om man har flere datapunkter än obekanta parametrar i funktionen.

Det finns en speciell typ av approximationsproblem där funktionen passerar genom datapunkterna. Detta innebär, att f (xi) = y_i, i = 1 : n. Vi säger i detta fall att f interpolerar data. Om approximations- funktionen är ett polynom kan vi uttrycka interpolationsproblemet p˚a följande sätt:

D˚a x₁, x₂, . . . , x_n och y₁, y₂, . . . , y_n är givna, sök ett polynom p_n−1(x) av gradtalet n − 1 (eller lägre) som uppfyller villkoret pn−1(x_i) = y_i för alla i = 1 : n.

S˚alunda interpolerar t.ex. polynomet p₂(x) = 1 + 4x − 2x² punkterna (−2, −15), (3, −5) och (1, 3) (se figuren).

(6)

Varje talpar (x_i, y_i) kan uppfattas som en punkt p˚a en kurva, beskriven av n˚agon funktion g(x) : y_i = g(x_i). Funktionen g kan vara explicit definierad, t. ex. om vi vill approximera sin x i punkterna x = 0, π/2 och π med ett andragradspolynom. Men funktionen kan ocks˚a vara implicit definierad, om vi t.ex. vill interpolera m¨atdata.

Hur uttrycks interpolanten p_n−1(x)? Istället för att använda de ”naturliga” baspolynomen 1, x och x², skulle vi kunna använda den alternativa basen 1, (x + 2) och (x + 2)(x − 3).

(7)

I denna bas kan det tidigare n¨amnda interpolationspolynomet p2(x) uttryckas p₂(x) = −15 + 2(x + 2) − 2(x + 2)(x − 3).

Olika baser har olika f¨ordelar.

När vi har bestämt oss för ett sätt att representera interpolanten, hur beräknar vi därp˚a koefficienterna?

Det visar sig att detta problem leder till l¨osningen av ett ekvationssystem med en koefficientmatris av en best¨amd symmetri.

Sedan koefficienterna blivit beräknade, hur kan vi p˚a ett effektivt sätt beräkna interpolanten? Detta problem klarnar när vi ritat en graf av polynomet.

Den klassiska metoden att ber¨akna ett interpolationspolynom har uppkallats efter Alexandre-Th´eophile Vandermonde (1735-1796, fransk matematiker, som bl.a. studerade determinanter). Som ett exempel skall vi studera problemet att finna det tredjegradspolynom

p₃(x) = a₁ + a₂x + a₃x² + a₄x³ som interpolerar punkterna (−2, 10), (−1, 4), (1, 6) och (2, 3).

(8)

Varje interpolationspunkt ger d˚a upphov till en line¨ar ekvation med fyra obekanta:

a₁ − 2a₂ + 4a₃ − 8a₄ = 10 a₁ − a₂ + a₃ − a₄ = 4 a₁ + a₂ + a₃ + a₄ = 6

a₁ + 2a₂ + 4a₃ + 8a₄ = 3 (1)

Vi kan uttrycka detta ekvationssystem i matrisform:







1 −2 4 −8

1 −1 1 −1

1 1 1 1

1 2 4 8











 a₁ a₂ a₃ a₄







=





 10

4 6 3







L¨osningen a=[4.5000 1.9167 0.5000 -0.9167]’ till dessa fyra ekvationer f˚ar man med f¨oljande enkla MATLAB-program:

y = [10; 4; 6; 3];

V = [1 -2 4 -8; 1 -1 1 -1; 1 1 1 1; 1 2 4 8];

a = V\y;

(9)

Vi skall nu studera det n–dimensionella interpolationsproblemet. I detta fall skall vi best¨amma koefficienterna a₁, a₂, . . . , a_n f¨or polynomet

p_n−1(x) = a₁ + a₂x + a₃x² + · · · + a_nxⁿ⁻¹ s˚a att

p_n−1(x_i) = a₁ + a₂x_i + a₃x²_i + · · · + a_nxⁿ⁻¹_i = y_i f¨or alla i = 1 : n. Genom att skriva ekvationerna i matrisform f˚ar vi







1 x₁ x²₁ · · · xⁿ⁻¹₁ 1 x₂ x²₂ · · · xⁿ⁻¹₂ 1 x₃ x²₃ · · · xⁿ⁻¹₃

... ... ... . . . ...

1 x_n x²_n · · · xⁿ⁻¹_n











 a₁ a₂ a₃ ...

a_n







=





 y₁ y₂ y₃ ...

y_n







Koefficientmatrisen (Vandermondes matris) kan vi kalla V . För att vi skall kunna lösa interpolationsproblemet, m˚aste V vara icke-singulär. Antag nu, att det finns en s˚adan vektor c, att V c = 0. Detta betyder, att polynomet

q(x) = c₁ + c₂x + · · · + c_nxⁿ⁻¹

(10)

skulle försvinna i var och en av de (disjunkta) punkterna x = x₁, x₂, . . . , x_n, m.a.o. q(x) är ett polynom av gradtalet n − 1 som har n skilda rötter. Detta kan endast inträffa om polynomet är ett noll–polynom (dvs c = 0). Allts˚a kan V inte vara en singulär matris.

Vi skall nu beskriva n˚agra sätt att konstruera matrisen V . Eftersom den i:te raden av V inneh˚aller potenser av x_i, och exponenterna växer fr˚an 0 till n − 1 d˚a man g˚ar fr˚an vänster till höger, s˚a kan man konstruera matriselementen med en vanlig dubbelslinga:

n = length(x); V = zeros(n,n);

for i = 1:n

% i:te raden:

for j = 1:n

V(i,j) = x(i)^(j-1);

end end

Detta ¨ar en radorienterad algoritm, eftersom den verkar p˚a matrisen rad f¨or rad.

Den inre slingan i denna skript kan vektoriseras med hj¨alp av elementvis exponentiering. S˚aledes ger t.ex.

MATLAB-kommandot u = [1 2 3 4].^[3 5 2 3] ˚at radvektorn u v¨ardet [1 32 9 64].

(11)

För att räkna ut elementen i den i:te raden av V m˚aste man upphöja skalären xi i var och en av exponenterna i radvektorn 0 : n − 1 = (0, 1, . . . , n − 1). Kommandot r = (x(i)*ones(1,n)).^(0:n-1) kommer därför att tillordna vektorn (1, x_i, x²_i, . . . , xⁿ⁻¹_i ) till r. Den i:te raden i V kan anges med V(i,:), och vi f˚ar d˚a skripten

for i=1:n

% i:te raden i V:

V(i,:) = (x(i)*ones(1,n)).^(0:n-1);

end

Vi kan ocks˚a kasta om slingorna i den ursprungliga skripten, varigenom vi f˚ar en kolumnorienterad algoritm:

for j=1:n

% j:te kolumnen:

for i = 1:n

V(i,j) = x(i)^(j-1);

end end

(12)

Om j > 1, s˚a kan V (i, j) uttryckas som produkten av x(i) och V (i, j − 1), vilket visar att potenseringarna kan g¨oras genom successiva multiplikationer:

V = ones(n,n);

for j=2:n

% j:te kolumnen:

for i = 1:n

V(i,j) = x(i)*V(i,j-1);

end end

Den j:te kolumnen bildas allts˚a genom elementvis vektormultiplikation:



 x₁

...

x_n



 . ∗





V_1,j−1 ...

V_n,j−1



 =



 V_1,j

...

V_n,j





som kan utföras med MATLAB–kommandot V(:,j) = x.*V(:,j-1). Vi kommer till sist fram till följande MATLAB-funktion för interpolering enligt Vandermondes metod:

(13)

function a = interpv(x,y)

%

% Invariabler:

% x: en kolumnvektor med n (olika) element

% y: en kolumnvektor med n element

%

% Utvariabel:

% a: en kolumnvektor med n element f¨or vilka g¨aller att om

% p(x) = a(1) + a(2)x + ... + a(n)x^(n-1) s˚a ¨ar

% p(x(i)) = y(i) , i = 1:n

%

n = length(x);

V = ones(n,n);

for j = 2:n

% j:te kolumnen:

V(:,j) = x.*V(:,j-1);

end

a = V\y;

(14)

Vi skall nu beskriva en effektiv metod för att beräkna polynomvärdena pn−1(x) = a₁ + . . . + a_nxⁿ⁻¹ för x = z, d˚a vi känner z och a(1:n). Vi kan beräkna värdet av pn−1(z) direkt genom att summera potenser av x:

n = length(a);

zk = 1;

pz = a(1);

for k = 2:n zk = z*zk;

pz = pz + a(k)*zk;

end

men det finns en effektivare metod. Vi skall för enkelhetens skull tillämpa den p˚a ett tredje gradens polynom, som vi skriver i formen p₃(x) = a₁ + a₂x + a₃x² + a₄x³ = ((a₄x + a₃)x + a₂)x + a₁. Polynomets värde i z kan följaktligen beräknas med programfragmentet

pz = a(4);

pz = z*pz + a(3);

pz = z*pz + a(2);

pz = z*pz + a(1);

(15)

F¨or ett godtyckligt gradtal kan man uttrycka denna kedjeregel med kommandosekvensen n = length(a);

pz = a(n);

for i=n-1:-1:1

pz = z*pz + a(i);

end

Denna metod, som kallas Horners regel ¨ar uppkallad efter William George Horner, en engelsk matematiker som beskrivit den 1819¹.

L˚at oss anta, att vi vill beräkna interpolanten i m˚anga, skilda punkter. Vi skall anta att vektorn z(1:m) har blivit initialiserad, och att vi vill tillordna värdet pn−1(z(i)) till pz(i) för alla i = 1 : m. Detta kan naturligtvis göras genom att man upprepar Horners metod i varje punkt:

1W.G. Horner,”A new method of solving numerical equations of all orders, by continuous approximation,”Phil. Trans., Vol. 109, 1819, ss. 308-335

(16)

m = length(z);

n = length(a);

pz = zeros(m,1);

for j = 1:m

% Ber¨akna p(z(j)).

pz(j) = a(n);

for i=n-1:-1:1

pz(j) = z(j)*pz(j) + a(i);

end end

Vi kan vektorisera detta programfragment genom att tänka p˚a vad det innebär att man ”samtidigt”beräknar interpolanterna i varje punkt z_i. Antag, att m = 5 och n = 4 (dvs vi vill beräkna en kubisk interpolant i fem punkter). Det första steget i Horners metod kan d˚a uttryckas







pz(1) pz(2) pz(3) pz(4)







=





 a(4) a(4) a(4) a(4)







(17)

I vektorform kan detta skrivas pz = a(n)*ones(m,1). Nästa steg kan göras simultant p˚a analogt sätt:







pz(1) pz(2) pz(3) pz(4)







=







z(1) ∗ pz(1) z(2) ∗ pz(2) z(3) ∗ pz(3) z(4) ∗ pz(4)





 +





 a(3) a(3) a(3) a(3)







eller i vektorform pz = z.*pz + a(3). F¨or en kubisk interpolant f˚ar vi allts˚a f¨oljande MATLAB-skript:

pz = a(4)*ones(m,1) pz = z.*pz + a(3);

pz = z.*pz + a(2);

pz = z.*pz + a(1);

Detta programfragment kan vi generalisera till följande MATLAB-rutin för att beräkna polynomvärden av vektorargument med Horners metod:

(18)

function pz = hornerv(a,z)

%

% Invariabler:

% a: en kolumnvektor med n element

% z: en kolumnvektor med m element

%

% Utvariabel:

% pz: en kolumnvektor med n element, f¨or vilka g¨aller att

% om p(x) = a(1) + ... + a(n)x^(n-1) , s˚a ¨ar

% pz(i) = p(z(i)) d˚a i=1:m.

%

n = length(a);

m = length(z);

pz = a(n)*ones(m,1);

for k = n-1:-1:1

pz = z.*pz + a(k);

end

Varje uppdatering av pz kr¨aver 2m flyttalsoperationer, s˚a att 2mn operationer beh¨ovs totalt.

(19)

Som en tillämpning skall vi studera ett program, som ritar kubiska interpolanter av sin x inom intervallet [0, 2π]. Punkternas x–koordinater väljs slumpmässigt med rand–funktionen, och sorteras i storleksordning med sort. Sinuskurvan anges med en heldragen linje, och interpolanten med en streckad linje.

% programfil: intsin

% Programmet ritar fyra slumpartade kubiska interpolanter f¨or

% sin(x) inom intervallet [0,2pi] och

% anv¨ander Vandermondes metod.

close

x1 = linspace(0,2*pi,100)’;

y1 = sin(x1);

for k=1:4

x = 2*pi*sort(rand(4,1));

y = sin(x);

a = interpv(x,y);

pz = hornerv(a,x1);

subplot(2,2,k)

plot(x1,y1,’-’,x1,pz,’--’,x,y,’o’) axis([0 2*pi -2 2])

end

(20)

Nedanst˚aende figur visar ett exempel p˚a grafiken, och visar ocks˚a hur valet av punkter p˚averkar ¨overensst¨am- melsen mellan den ursprungliga kurvan och interpolanten.