Från förra gången: Newton-Raphsons metod
Idé:
För att hitta en rot till f(x)=0
utgår man från en första
Approximation x0 och
använder derivatan för att
dra en tangent som skär
x-axeln närmare roten och
upprepar detta tills man är
tillräckligt nära:
xn+1= xn – f(xn)/f’(xn)
Om f är ett formeluttryck
så kan man normalt lätt
bilda derivatan f’
Yngve Sundblad
Föreläsning 7 sid.1
SF 1518/19 ht 2015
ACPU - 29 January 2004
21 sept.
Yngve Sundblad
Exempel (EXS 4.6) med Newton-Raphson
P(x)=4x4 - 7x3 – 5.5x2 +27.5x – 50 = 0. Hitta noga roten nära 1.5.
P’(x)=16x3 – 21x2 – 11x +27.5 p=[4 -7 5.5 27.5 -50]
for k=1:4 %derivatan
pprim(k)=(5-k)*p(k);
end
x=1.5; dx=1;
format long
while abs(dx/x)>1e-12
px=polyval(p,x);
pprimx=polyval(pprim,x);
dx=-px/pprimx;
disp(x); disp(dx);
x=x+dx;
end
ger
1.500000000000000 -0.004926108374384
1.495073891625616 -0.000013421067689
1.495060470557926 -0.000000000099026
1.495060470458900 0.000000000000000
Yngve Sundblad
Föreläsning 7 sid.2
Fördubblat antal
decimaler varje varv Kvadratisk konvergens?
SF 1518/19 ht 2015
ACPU - 29 January 2004
21 sept.
Yngve Sundblad
Nytt: Trunkationsfel för Newton-Raphson
Låt roten till f(x)=0 vara r.
Taylorutveckla:
0=f(r)=f(xn +(r-xn))=f(xn)+(r-xn)*f’(xn)+(r-xn)2*f’’(xn)/2 + …
Medelvärdessatsen: Om vi avbryter ersätts f’’(xn) med f’’(µ),
där µ ligger mellan xn och r.
Dividera med f’(xn) och stuva om xn – f(xn)/f’(xn) – r = (r-xn)2*f’’(µ)/(2*f’(xn))
xn+1–r = (r-xn)2*f’’(µ)/2/f’(xn) eller när n går mot ∞ konvergerar
(xn+1–r)/(r-xn)2 mot K = f’’(r)/(2*f’(r))
Alltså kvadratisk konvergens!
Yngve Sundblad
Föreläsning 7 sid.3
SF 1518/19 ht 2015
ACPU - 29 January 2004
21 sept.
Yngve Sundblad
Att hitta startvärde för ekvationslösning
Ofta ger en graf tillräcklig information men den kan luras. F=x–4sin(2x)–243/80=0 plottas [–1:8] (rot>–4+3.03, rot<4+3.04)
x=-1:0.02:8;
Förstorad
f=x-4*sin(2*x)-243/80;
plot([-1 8],[0 0],x,f)
x=6.98:0.0001:7.02;
f=x-4*sin(2*x)-243/80; plot([6.98 7.02],[0 0],x,f)
Förstoringen lurades pga. Begränsad linjelängd! EXS3.7 liknande, analys krävs!
Yngve Sundblad
Föreläsning 7 sid.4
SF 1518/19 ht 2015
ACPU - 29 January 2004
21 sept.
Yngve Sundblad
Sekantmetoden
När f(x)är svår eller omöjlig att derivera är
sekantmetoden ett bra alternativ.
Man drar alltså en sekant ,mellan punkten (a,f(a)) och (b,f(b)) och
låter skärningspunkten vara nya (bättre) x och upprepar detta tills intervallet är tillräckligt
litet.
Nya punkten xn+2 beräknas ur de två tidigare med
xn+2 =xn+1–(xn+1–xn)/(f(xn+1)–f(xn))*f(xn+1)
Yngve Sundblad
Föreläsning 7 sid.5
SF 1518/19 ht 2015
ACPU - 29 January 2004
21 sept.
Yngve Sundblad
Exempel (EXS 4.6) med sekantmetoden
P(x)=4x4 - 7x3 – 5.5x2 +27.5x – 50 = 0. Hitta noga roten nära 1.5.
Format long; x0=1; x1=2;it=1;
P=[4 -7 5.5 27.5 -50];
Ger
0.636363636363636
-0.095652128035623
-0.038316357619555
0.002594303021054
-0.000049854344873
-0.000000069845456
0.000000000001916
Iterationer:9 Rot: 1.495060470459
f0=polyval(P,x0); h=1;
while abs(h)/abs(x1)>1e-12
f1=polyval(P,x1);
h=(x1-x0)/(f1-f0)*f1;
disp(h);
x0=x1;f0=f1;
x1=x1-h; it=it+1;
end
fprintf('Iterationer: %d \n', it);
fprintf('Rot: %14.12f \n', x1)
Yngve Sundblad
Föreläsning 7 sid.6
Knapp fördubbling av decimalerna. Man kan visa att osäkerheten ~K*h1.6.
SF 1518/19 ht 2015
ACPU - 29 January 2004
21 sept.
Yngve Sundblad
Fixpunktiteration
Skriv om f(x)=0 på formen x=G(x). Vi letar efter ”fixpunkten” där x nära avbildas på sig själv. Iterera xn+1=G(xn) med x0 som startvärde, skattning av roten.
Om r är roten så xn+1–r = G(xn) – G(r) = G’(µ)*(xn–r) (medelvärdessatsen)
(xn+1–r )/(xn+1–r) = G’(µ) ≈ m
Om kring roten gäller |m| < 1 så konvegerar den linjärt med m
gånger mindre osäkerhet för varje iteration.
Exempel (EXS 2.14): En 400 meters ellipsformad löparbana ska
anläggas på en 160 m lång plan. Hur bred plan krävs?
För omkretsen av en ellips använder vi en formel av självlärda
indiska matematiksnillet Srinivasa Ramanujan (1887-1920): π(a+b)(1+3c/(10√(4-3c)), där c=(a-b)2/(a+b)2; Fotnot: a=b ger 2πa, cirkeln
Yngve Sundblad
Föreläsning 7 sid.7
SF 1518/19 ht 2015
ACPU - 29 January 2004
21 sept.
Yngve Sundblad
Löparbanan
Vi söker alltså b så att π(a+b)(1+3c/(10√(4-3c))=400, där c=(a-b)2/(a+b)2, a=80
Skriv om för iteration: b=400/(π*(1+3c/10√(4-3c)))–80
Startvärde (rektangel): b=40.
function res=biterat(a,b)
format short;
b0=40; it=1; b1=biterat(80,b0)
while abs(b1-b0)>0.05
it=it+1;
b0=b1;
b1=biterat(80,b0)
end
fprintf('Bredd: %f \n',b1)
fprintf('efter %d iterationer \n', it)
c=(a-b)^2/(a+b)^2;
res=400/pi/(1+3*c/
(10+sqrt(4-3*c)))-a;
ger 43.8588 44.6561 44.8029 44.8292 Bredd:44.83 efter 4 iterationer
end
m=G’(r) ≈ 0.2 (0.03/0.15)
Yngve Sundblad
Föreläsning 7 sid.8
SF 1518/19 ht 2015
ACPU - 29 January 2004
21 sept.
Yngve Sundblad
Icke-linjära ekvationssystem
Modifierade varianter av metoderna för att lösa ekvationer kan
användas för att lösa system av sådana. Vi tar upp
fixpunktiteration och Newton-Raphson.
Allmänt kan ett ekvationssystem med n ekvationer i n variabler skrivas f1(x1, x2, x3, … , xn)=0 f2(x1, x2, x3, … , xn)=0 … fn(x1, x2, x3, … , xn)=0
Vänsterledets ”derivata” är Jacobianmatrisen: ∂f1/∂x1 ∂f1/∂x2 …. ∂f1/∂xn
J(x) = (df/dx),
∂f2/∂x1 ∂f2/∂x2 …. ∂f2/∂xn
specialfall: J=
…
Ax–b=0, J=A
∂fn/∂x1 ∂fn/∂x2 …. ∂fn/∂xn
f(x)=0, J=df/dx
(
Yngve Sundblad
)
Föreläsning 7 sid.9
SF 1518/19 ht 2015
ACPU - 29 January 2004
21 sept.
Yngve Sundblad
Newton-Raphson för ekvationssystem
För att hitta en rot till f(x)=0 använder man motsvarande formel
som i envariabel- (skalära) fallet: Beräkna successiva c(n) ur ekvationssystemet J(x(n))*c(n) = f(x(n))
Sätt x(n+1) = x(n) – c(n)
Avbryt när || c(n) || < eps, feltoleransen (tillåten osäkerhet)
( Jämför med skalära xn+1= xn – f(xn)/f’(xn) )
Konvergensen blir som i skalära fallet kvadratisk || h(n+1) || / || h(n) ||2 går mot K när n går mot ∞.
Beviset är analogt med skalära fallet.
Yngve Sundblad
Föreläsning 7 sid.10
SF 1518/19 ht 2015
ACPU - 29 January 2004
21 sept.
Yngve Sundblad
Exempel (≈ EXS3.9)
Vi vill hitta skärningspunkterna mellan lemniskatan (x2 + y2)2 = x2 – y2 och parabeln y = x2 – 1/2
Vi byter till polära koordinater: x=r*cosv, y=r*sinv, sätt in: r4=r2 (cos2v–sin2v), r = √cos(2v), rita med parabeln: v=[0:pi/100:2*pi]; r=sqrt(cos(2*v)); polar(v,r);hold on; x=-1:0.01:1;f=x.2-1/2; plot([-1 1],[0 0],x,f)
Skärningspunkter ≈ (±0.8,0.2), (±0.4,–0.35)
Warning: Imaginary parts of complex X and/or Y arguments ignored > In polar (line 192)
Gick bra ändå!
Yngve Sundblad
Föreläsning 7 sid.11
SF 1518/19 ht 2015
ACPU - 29 January 2004
21 sept.
Yngve Sundblad
Newton-Raphson på lemniskata/parabel
(x2 + y2)2 – x2 + y2 = 0,
Jacobianen: J =
(
y – x2 + ½ = 0
2x*2(x2 + y2)–2x
–2x
)
2y*2(x2 + y2)+2y 1
x0=[0.4 0.8]; y0=[-0.35 0.2]; dcnorm=1; varv=0;
for k=1:2
x=x0(k); y=y0(k); c=[x y]’; while abs(dcnorm)>5e-4
f=[(x^2+y^2)^2-x^2+y^2 y-x^2+1/2]'
J=[4*x*(x^2+y^2)-2*x 4*y*(x^2+y^2)+2*y;
-2*x 1]
ger (med beaktande av ±x): dc=-J\f; c=c+dc
c = [x y] = [±0.4260 -0.3185] dcnorm=norm(dc,Inf)
x=c(1); y=c(2); varv=varv+1; dcnorm = 1e-6 c = [±0.8730 0.2636]
end dcnorm = 8e-6
end
Yngve Sundblad
Föreläsning 7 sid.12
SF 1518/19 ht 2015
ACPU - 29 January 2004
21 sept.
Yngve Sundblad
Fixpunktsmetoder för ekvationssystem
Om man skriver om ekvationssystemet som x = G(x) och
jakobianmatrisens norm, M=||dG/dx||, är mindre än 1, ju mindre dess
bättre, kan man använda Jacobis eller Gauss-Seidels metod, som vi
gick igenom och använde på linjära system.
Anledningen är ofta att det är svårt att uppskatta jakobianen och man
använder uppskattningen M≈||x(n+1)–x(n)||/||x(n)–x(n-1)||
Felskattningen (osäkerheten) lösningen r kan härledas till ||x(n+1)–r|| ≈ M*||x(n+1)–x(n)||/(1-M)
Om man kan skriva om ekvationssystemet med en linjär del: f(x)= Ax + g(x) = 0, så kan man ofta använda Picarditeration:
A x(n+1) = – g(x(n)) med ett bra startvärde, x(0), ger lösningen med
tillräcklig noggrannhet, med trunkationsfelet uppkattat som ovan.
Yngve Sundblad
Föreläsning 7 sid.13
SF 1518/19 ht 2015
ACPU - 29 January 2004
21 sept.
Yngve Sundblad
Exempel på Picarditeration
Vi har ekvationssystemet (omskrivet med linjära termerna till vänster): 10x1 – x2 = x12 + x22 x1 + 10x2 – x3 = 2 – x23; x1 + 3 x3 = 1 – x33
Ax = g(x), A = [10 –1 0; 1 10 –1; 1 0 3] g(x)=[x12+x22; 2–x2; 1–x33]
Eftersom vänsterledet är så diagonaltungt är värdena i lösningen små. För att få startvärde ignorerar vi de icke-linjära termerna och löser Ax(0) = [0 2 1]’ Programmet ger x0=[0.023 0.230 0.326], x=[0.027989 0.227400 0.313714] trunk=1.5e-6 varv=4
Yngve Sundblad
A=[10 -1 0; 1 10 -1; 1 0 3]; b=[0 2 1]';
x=A\b
dx=[1 1 1]; varv=0; trunk=1;
format long;
while trunk>1e-7
xold=x;olddx=dx;
g=[x(1)^2+x(2)^2 2-x(2)^3 1-x(3)^3]’;
x=A\g;
dx=x-xold;
M=norm(dx)/norm(olddx);
trunk=M*norm(dx)/(1-M)
varv=varv+1;
end
disp(x); disp(trunk); disp(varv)
Föreläsning 7 sid.14
SF 1518/19 ht 2015
ACPU - 29 January 2004
21 sept.
Yngve Sundblad
Ekvationer och –system i kurslitteraturen
Det vi gått igenom (och ingår i kursen) svarar mot
stoffet i Pohl i hela kapitel 3 och 4 utom avsnittet
4:2C Ickelinjära minstakvadratproblem.
I NAM (Eriksson) motsvarar det vi gått igenom hela stoffet i
kapitel 6.
Yngve Sundblad
Föreläsning 7 sid.15
SF 1518/19 ht 2015
ACPU - 29 January 2004
21 sept.
Yngve Sundblad