1. No category

Institutionen för datavetenskap
Umeå universitet
25 mars 2015
Teknisk beräkningsvetenskap I 5DV154
Deltentamen inkusive svar
Tid: 15.00–19.00
Hjälpmedel: Matlab.
Maximalt antal poäng: 100
50 poäng är tillräckligt för godkänt på detta delprov
1.
(a) När man löser ekvationssystemet Ax = b i Matlab, så är det rekommenderat att skriva
x=A\b och inte x=inv(A)*b. Varför?
[2p]
(b) Nämn en fördel och en nackdel med att använda Newtons metod istället för fixpunktsiteration för att lösa ickelinjära ekvationer.
[2p]
(c) Låt (x 0 , y 0 ), . . . , (x n , y n ) vara en mängd av talpar. (i) Ge ett villkor för att punktmängden
ska kunna interpoleras med polynom. (ii) Vilka problem kan uppstå när man interpolerar
godtyckliga punktmängder med polynom?
[3p]
(d) Ge villkor för att den skalära, linjära ordinära differential ekvationen y = λy + f (t ) ska
vara stabil, asymptotiskt stabil, samt ostabil med avseende på begynnelsevärden.
[3p]
0
2.
(a) Skriv en Matlabfunktion som löser ekvationssystemet f i (x) = 0, i = 1, 2, . . . , N med
Newtons metod. Funktionshuvudet är givet nedan och funktionen ska fungera enligt
den specifikation som ges i detta.
function root = newton(fun,x0,tol)
%NEWTON numerically solves equation system by using Newton’s mehod
%
%
ROOT = NEWTON(FUN,X0,TOL) find an approximate solution ROOT to
%
the equation system FUN(X)=0 f by using Newton’s method with XO
%
as starting point and TOL as the stopping tolerance (that is,
%
NORM(FUN(ROOT)) <= TOL). FUN is a function handle. The function
%
[Y,J]=FUN(X) takes a vector argument X and returns the function
%
value Y and the Jacobian J.
[12p]
(b) Vi studerar en pendel som som svänger med amplitud 1 och är intresserade av att
bestämma frekvens och fas. Positionen av pendeln vid tiden t ges av y(t ) = sin(ωt + ϕ),
där ω och ϕ är reella konstanter. Antag att vi har följande mättdata
t 1 = 0.2, y 1 = 0.65, t 2 = 0.4, och y 2 = 0.90.
(i) Sätt upp problemet att bestämma de okända konstanterna a 1 och a 2 som ett system
på formen f(a) = 0, där a = (ω, ϕ)T, samt beräkna Jakobianmatrisen J av f.
[6p]
(ii) Skriv en Matlab function pendel som tar a = (ω, ϕ)T som inargument och beräknar och returnerar funktionsvärdena f samt Jacobianmatrisen J för det problem vi
studerar i denna uppgift.
[6p]
(c) Hitta konstanterna a genom att använda dig av dina Matlab funktioner newton och
pendel samt startgissningen a0 = [0, 1]T.
[6p]
3.
(a) Givet den generella kvadraturformeln
Z
1
−1
f (x) dx ≈ α1 f (−1) + α2 f (0) + α3 f (1),
härled villkor på konstanterna α1 , α2 och α3 så att formeln är exakt för alla
(i) konstanta funktioner f ,
(ii) linjära polynom f ,
(iii) kvadratiska polynom f .
Vad heter denna metoden som du härlett?
[10p]
(b) Beskriv hur man från den metod som du härlett
R 1 ovan kan skapa en så kallad sammansatt
kvadraturmetod för att beräkna integralen 0 f (x) dx.
R1
(c) Skriv en Matlabfunktion som beräknar integralen 0 f dx med den sammansatta metoden du härlett i (a–b). Funktionshuvudet är givet nedan och funktionen ska fungera
enligt den specifikation som ges i detta.
[10p]
function I = myquad(fun,M)
%MYQUAD approximates integral by using a composite quadrature rule
%
%
I = MYQUAD(FUN,M) approximates the integral of scalar-valued
%
function FUN from 0 to 1 by an appropriate quadrature rule that
%
evaluates the function FUN in M points. FUN is a function handle
%
and M is an odd positive integer. The function Y=FUN(X) should
%
accept a vector argument X and return a vector result Y, the
%
integrand evaluated at each element of X.
(Om du inte löst deluppgift (a–b) så använd en av de sammansatta kvadraturmetoderna
som ingår i kursen.)
[10p]
(d) Bestäm numeriskt metoden noggrannetsordning då den används för att beräkna
1
Z
I=
ex dx
0
(Observera att vi kan räkna ut integralen ovan exakt: I = e − 1.)
4.
[5p]
(a) Begynnelsevärdesproblemet för y 0 = f (t , y) kan lösas numeriskt med trapetsmetoden
y k+1 = y k +
¢
∆t ¡
f (t k+1 , y k+1 ) + f (t k , y k ) .
2
(i) Är denna metod implicit eller explicit?
[2p]
(ii) Bestäm metodens stabilitetsvillkor på negativa reella axeln, d.v.s. villkoret för stabilitet när schemat tillämpas på problemet y 0 = λy för λ < 0.
[13p]
(b) Sätt upp följande begynnelsevärdesproblem på standardform
y 00 = −t y − sin(y 0 )| sin(y 0 )|
y(0) = 1, y 0 (0) = 1
och lös detta problem från tiden 0 till tiden T = 5 med med Matlabs inbyggda ode lösare
ode45 samt plotta lösningen y som funktion av t .
[10p]
Svar 1.
(a) När man skriver A\b så löses det linjära systemet med Gausselimination, vilket tar
mycket mindre antal flyttalsoperationer jämfört med att explicit beräkna inversen av A
och därefter utföra matris–vektormultiplikationen inv(A)*b.
(b) Fördel: snabbare konvergens, nackdel: behöver derivator
(c) (i) Om alla xi är skilda kan punkterna alltid interpoleras. (ii) Interpolation med ett
polynom av hög grad kan medföra stora oscillationer mellan interpolationspunkterna.
(d) En skalär, linjär ODE skriven på formen y 0 = λy + f (t ) är stabil när ℜ(λ) ≤ 0, asymptotisk
stabil när ℜ(λ) < 0 och ostabil annars (alltså om ℜ(λ) > 0). Således är (a) och (b) instabila,
(c) är asymptotiskt stabil och (d) är stabil.
Svar 2.
(a) function root = newton(fun,x0,tol)
%NEWTON numerically solves equation system by using Newton’s mehod
%
%
ROOT = NEWTON(FUN,X0,TOL) find an approximate solution ROOT to
%
the equation system FUN(X)=0 f by using Newton’s method with XO
%
as starting point and TOL as the stopping tolerance (that is,
%
NORM(FUN(ROOT)) <= TOL). FUN is a function handle. The function
%
[Y,J]=FUN(X) takes a vector argument X and returns the function
%
value Y and the Jacobian J.
root = x0;
[fval,J] = fun(root);
iter = 0;
maxiter = 100;
while norm(fval)>tol && iter<maxiter
s = -J\fval;
root = root + s;
[fval,J] = fun(root);
iter = iter+1;
end
Eller alternativt (om man inte vill begränsa antalet iterationer...)
function root = newton(fun,x0,tol)
%
[ ... ]
root = x0;
[fval,J] = fun(root);
while abs(fval)>tol
s = -J\fval;
root = root + s;
[fval,J] = fun(root);
end
(b) Problemet som vi vill lösa är: Finn ω och ϕ så att
sin(ωt 1 + ϕ) − y 1 = 0,
sin(ωt 2 + ϕ) − y 2 = 0.
Ovanstående system kan skrivas på formen f(a) = 0, där a = (ω, ϕ)T och
µ
¶
sin(a 1 t 1 + a 2 ) − y 1
f(a) =
.
sin(a 1 t 2 + a 2 ) − y 2
Jacobianen av f ges av
 ∂f
J=
∂a 1
∂ f1
∂a 2
∂ f2
∂a 1
∂ f2
∂a 2
1

=
Ã
t 1 cos(a 1 t 1 + a 2 )
cos(a 1 t 1 + a 2 )
t 2 cos(a 1 t 2 + a 2 )
cos(a 1 t 2 + a 2 )
function [f,J] = pendel(a)
t = [0.2; 0.4];
y = [0.65; 0.90];
f = [sin(a(1)*t(1)+a(2)) - y(1);
sin(a(1)*t(2)+a(2)) - y(2)];
J = [t(1)*cos(a(1)*t(1)+a(2)) cos(a(1)*t(1)+a(2));
t(2)*cos(a(1)*t(2)+a(2)) cos(a(1)*t(2)+a(2))];
(c) En lösning fås genom att göra anropet
>> root = newton(@pendel,[0; 1],1e-5)
root =
2.0609
0.2954
!
.
Svar 3.
(a) (i) Ansätts f (x) = a får vi
1
Z
−1
a dx = 2a = a(α1 + α2 + α3 ),
vilket ger villkoret
α1 + α2 + α3 = 2
(1)
för att formeln skall vara exakt för alla konstanta funktioner. (ii) Ansätts f (x) = bx får vi
(obs att f är en udda funktion!)
Z
1
−1
bx dx = 0 = b(−α1 + α3 ),
vilket ger villkoret
α1 = α3 .
(2)
Formeln är exakt för alla linjära funktioner f (x) = a + bx om villkoren (1) och (2) bägge
är uppfyllda. (iii) Ansätts f (x) = cx 2 får vi
Z
1
−1
cx 2 dx = c
2
= c(α1 + α3 ),
3
d.v.s.
α1 + α3 =
2
3
(3)
Formeln är exakt för alla kvadratiska funktioner f (x) = a + bx + cx 2 om villkoren (1), (2)
och (3) alla är uppfyllda. Tillsammans ger dessa villkor att
1
α1 = α3 = ,
3
4
α2 = ,
3
vilket är Simpsons formel.
(b) För att göra en sammansatt kvadraturregel, delar vi upp intervallet som vi vill integrera i
delintervall och använder oss av en skalad “enkel” kvadraturregel på varje intervall. Det
första steget är att dela upp intervallet [0, 1] i N delintervall [x i , x i +1 ], i = 0, .R. . , N − 1, där
1
x i = i ∆x, i = 0, . . . , N och ∆x = 1/N . Från kvadraturformeln ovan har vi att −1 f (x) dx ≈
( f (−1) + 4 f (0) + f (1))/3; vilket medför att för varje delintervall så har vi att
Z
x i +1
xi
f (x) dx ≈ ∆x
f (x i ) + 4 f ((x i + x i +1 )/2) + f (x i +1 )
6
Genom att summera över alla delintervall får vi att
1
Z
0
f (x) dx =
NX
−1 Z x i +1
i =0
xi
f (x) dx ≈
NX
−1
i =0
∆x
f (x i ) + 4 f ((x i + x i +1 )/2) + f (x i +1 )
6
En alternativ lösning använder sig av tittar istället på M = 2N + 1 punkter x̂ i = i ∆x̂, där ,
ˆ = 1/(2N ), vilket ger slututtrycket:
i = 0, . . . , M − 1 och ∆x
1
Z
0
f (x) dx ≈
´
∆x̂ ³
f (x̂ 0 ) + 4 f (x̂ 1 ) + 2 f (x̂ 2 ) + 4 f (x̂ 3 ) + . . . + 2 f (x̂ M −2 ) + 4 f (x̂ M −1 ) + f (x̂ M )
3
(c) Observera att här ska vi använda oss av M punkter, x och dx i koden nedan motsvarar x̂
och ∆x̂ i lösningen ovan
function I = myquad(fun,M)
%MYQUAD approximates integral by using a composite quadrature rule
%
%
I = MYQUAD(FUN,M) approximates the integral of scalar-valued
%
function FUN from 0 to 1 by an appropriate quadrature rule that
%
evaluates the function FUN in M points. FUN is a function handle
%
and M is an odd positive integer. The function Y=FUN(X) should
%
accept a vector argument X and return a vector result Y, the
%
integrand evaluated at each element of X.
x = linspace(0,1,M);
dx = x(2);
weights = ones(size(x));
weights(2:2:end-1) = 4;
weights(3:2:end-2) = 2;
weights = weights*dx/3;
fval = fun(x);
I = fval*weights’;
(d) Låt
>> qfun = @(x) exp(x);
Testa först med ∆x = 1/2, 1/4, 1/8 samt 1/16 vilket ger 3, 5, 9 respektive 17 punkter. Felet
i integralberäkningen ges i dessa fall av
>> error3 = myquad(qfun,3)-exp(1)+1
error3 =
5.7932e-04
>> error5 = myquad(qfun,5)-exp(1)+1
error5 =
3.7013e-05
>> error9 = myquad(qfun,9)-exp(1)+1
error9 =
2.3262e-06
>> error17 = myquad(qfun,17)-exp(1)+1
error17 =
1.4559e-07
Genom att beräkna
>> error3/error5
ans =
15.6517
>> error5/error9
ans =
15.9113
>> error9/error17
ans =
15.9777
ser vi att felet minskar med en faktor ungefär 16 när vi halverar ∆x vilket tyder på att
noggrannhetsordningen är 4 (eftersom 24 = 16).
Svar 4.
(a)
(i) Metoden är implicit
(ii) För y 0 = λy ger schemat att
y k+1 = y k +
¢
∆t ¡
λy k+1 + λy k
2
vilket ger oss
(1 − ∆t λ/2)y k+1 = (1 + ∆t λ/2)y k
För stabilitet behöver vi |y k+1 | ≤ |y k |, vilket gäller om och endast om
¯
¯
¯1 + ∆t λ/2¯
≤ 1.
|1 − ∆t λ/2|
Eftersom λ < 0, kan vi skriva λ = −|λ| och multiplicera båda sidor av olikheten ovan
med |1 − ∆t λ/2| = 1 + ∆t |λ| ≥ 0 för att erhålla
¯
¯
¯1 − ∆t |λ|/2¯ ≤ 1 + ∆t |λ|/2,
vilket ger följande två olikheter
−1 − ∆t |λ|/2 ≤ 1 − ∆t |λ|/2 ≤ 1 + ∆t |λ|/2.
Dessa olikheter är alltid uppfyllda (ty ∆t |λ| ≥ 0) och således är schemat ovillkorligt
stabilt.
(b) Låt u 1 = y och u 2 = y 0 . Vi har då att
u0 = f(t , u) =
µ
u2
−t u 1 − sin(u 2 )| sin(u 2 )|
¶
u(0) =
µ ¶
1
1
>> testode = @(t,u) [u(2); -t*u(1)-sin(u(2))*abs(sin(u(2)))];
>> [t,u] = ode45(testode,[0 5],[1; 1]);
>> plot(t,u(:,1))