Projekt: Inovace výuky optiky se zaměřením na získání experimentálních dovedností Registrační číslo: CZ.1.07/2.2.00/28.0157
Numerické metody a programování Lekce 1
Tento projekt je spolufinancován Evropským sociálním fondem a státním rozpočtem České republiky.
Numerické metody a programování
Obsah přednášky 1. Mathematica: základy programování, symbolické výpočty, vizualizace dat. 2. Programování v prostředích Matlab/Octave: srovnání s jazykem C, knihovny funkcí, vytváření uživatelských funkcí. 3. Úvod do numerických metod: přesnost, zaokrouhlovací chyby, stabilita. 4. Algebra: práce s vektory a maticemi, řešení algebraických rovnic, SVD, Choleského dekompozice. 5. Aproximace funkcí: interpolace, extrapolace, interpolační mnohočleny, spliny. 6. Numerická integrace/derivace: elementární a pokročilé algoritmy, multidimenzionální integrace, integrace obyčejných diferenciálních rovnic. 7. Řešení soustav nelineárních rovnic: bisekce, NewtonovaRaphsonova metoda. 8. Optimalizace: gradientní metody, “downhill simplex” ve více dimenzích, metoda konjugovaného gradientu, lineární programování. 9. Modelování: metoda nejmenších čtverců, teorie odhadu, nelineární modely, konfidenční intervaly. 10. Fourierova transformace: spojitá a diskrétní transformace, FFT algoritmus a jeho použití, Nyquistova frekvence, diskrétní Fourierova transformace 2D a 3D. 11. Aplikace I: numerická simulace šíření optického signálu, Fresnelova difrakce, požadavky na vzorkování, aliasing. 12. Aplikace II: analýza zobrazovacích systémů a aberací, šíření v turbulentním prostředí, rekonstrukce vlnoplochy.
Doporučená literatura • • • •
E. Vitásek, Numerické metody (SNTL, Praha, 1987). W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes in C, (Cambridge University Press, Cambridge, 1992); dostupné online na http://www.nr.com J.D. Schmidt, Numerical Simulation of Optical Wave Propagation (SPIE Press, 2010) Manuály Matlab/Octave (http://www.octave.cz), Mathematica, Oslo (http://www.lambdares.com)
Mathematica (Wolfram Research) http://www.wolfram.com/mathematica/ • • •
symbolické výpočty numerické výpočty vizualizace dat a výsledků
H* prirazeni, relace *L a=1 1 a 1 a =. a a a=1 1 a = a+1 2 a ++; a 3 a += 5 8 b=1 1 a>b True ab False
H* komplexni cisla *L x = 2+I 2+ä x^2 3+4ä Re@xD 2
1
math_tisk.nb
Im@xD 1 Abs@xD 5 Arg@xD 1 ArcTanA E 2 Conjugate@xD 2-ä x =. f = x ^ 2 - Abs@xD ^ 2 x2 - Abs@xD2 Simplify@fD x2 - Abs@xD2 Simplify@f, Im@xD == 0D 0
H* ridici struktury *L a = 2; b = 4; If@a < b, mensi = a, mensi = bD; mensi 2 For@i = 1, i £ 10, i ++, Print@iDD 1 2 3 4 5 6 7 8 9 10
2
math_tisk.nb
suma = 0 0 For@i = 1, i £ 10, i ++, suma += iD Print@sumaD 55 suma = 0; i = 1; While@i £ 10, suma += i; i += 2D Print@sumaD 25 min@a_, b_D := If@a < b, a, bD min@1, 2D 1 min@3, 4D 3
H* vektory, matice *L v = 81, 2, 3< 81, 2, 3< Sqrt@vD 2 ,
91,
3 =
v.v 14 v = Table@Cos@xD, 8x, 0, 2 Π, Π 2
Π , Π,
90, 2
, 2 Π= 2
Cos@v2D 81, 0, - 1, 0, 1< m = 882, 1<, 8- 1, 3<< 882, 1<, 8- 1, 3<< MatrixForm@mD K
2 1 O -1 3
3
math_tisk.nb
m@@1, 2DD 1 v = 881<, 81<< 881<, 81<< MatrixForm@vD K
1 O 1
[email protected] K
3 O 2
a = 88a11, a12<, 8a21, a22<< 88a11, a12<, 8a21, a22<< b = 88b11, b12<, 8b21, b22<< 88b11, b12<, 8b21, b22<< c = a.b; c MatrixForm K
a11 b11 + a12 b21 a11 b12 + a12 b22 O a21 b11 + a22 b21 a21 b12 + a22 b22
Det@aD - a12 a21 + a11 a22 Tr@aD a11 + a22
H* linearni algebra *L a = 880, 1, 1<, 81, 1, 1<, 81, 1, 0<<; a MatrixForm 0 1 1 1 1 1 1 1 0 eig = Eigenvalues@aD 91 +
2 , - 1, 1 -
2 =
vec = Eigenvectors@aD 991,
2 , 1=, 8- 1, 0, 1<, 91, -
2 , 1==
a.vec@@1DD 91 +
2 , 2+
2 , 1+
2 =
4
math_tisk.nb
eig@@1DD vec@@1DD 91 +
2 ,
2 J1 +
2 N, 1 +
2 =
b = Inverse@aD; b MatrixForm -1 1 0 1 -1 1 0 1 -1 a.b MatrixForm 1 0 0 0 1 0 0 0 1 prava = 81, 2, 3< 81, 2, 3< x = b.prava 81, 2, - 1< a.x prava True
H* vyrazy *L a =.; b =.; c =.; x =.; f1 = x x f2 = Exp@- xD ã-x f = f1 * f2 ã-x x
H* derivace *L D@f, xD ã-x - ã-x x D@f, x, xD - 2 ã-x + ã-x x
H* integrace *L
5
math_tisk.nb
integral = Integrate@f, xD ã-x H- 1 - xL tem = D@integral, xD - ã-x - ã-x H- 1 - xL simp = Simplify@temD ã-x x simp == f True Integrate@Exp@- x ^ 2D, 8x, - Infinity, Infinity
2
, IntegrateAã-a x , 8x, - ¥, ¥<, Assumptions ® Re@aD £ 0EE
IfARe@aD > 0, a
Simplify@vysl, a > 0D Π a Plot@Cos@x ^ 2D, 8x, 0, 5
0.5
1
2
3
4
5
-0.5
-1.0
Integrate@Cos@x ^ 2D, 8x, 0, Α
Π FresnelCA 2
ΑE Π
6
math_tisk.nb
Plot@Integrate@Cos@x ^ 2D, 8x, 0, Α
0.8
0.6
0.4
0.2
0
2
4
6
8
10
Integrate@Cos@x ^ 2D, 8x, 0, 1
NIntegrate::ncvb : NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in x near 8x< = 80.<. NIntegrate obtained 0.4808197282114536` and 0.0116887984970139` for the integral and error estimates. 0.48082 i2 = NIntegrate@Cos@x ^ 2D, 8x, 0, 100<, MaxRecursion ® 10D 0.625129 i1 - i2 8.54872 ´ 10-15
H* soucty rad *L Sum@n ^ 2, 8n, 1, 3
7
math_tisk.nb
Sum@1 2 ^ n, 8n, 0, Infinity
H* rovnice *L Solve@2 x + 5 0, xD 5 99x ® -
== 2
Solve@82 x + y 1, x - y == 2<, 8x, y
b2 - 4 a c
99x ®
-b +
b2 - 4 a c
=, 9x ® 2a
== 2a
f = Expand@Hx - 2L Hx - 1L Hx + 2LD 4 - 4 x - x2 + x3 Solve@f 0, xD 88x ® - 2<, 8x ® 1<, 8x ® 2<< f = Cos@xD - x - x + Cos@xD Solve@f 0, xD
Solve::tdep : The equations appear to involve the variables to be solved for in an essentially non-algebraic way. Solve@- x + Cos@xD 0, xD vysl = FindRoot@f, 8x, - 1
H* diferencialni rovnice *L
8
math_tisk.nb
DSolve@y '@xD + y@xD a Sin@xD, y@xD, xD 1 99y@xD ® ã-x C@1D +
a H- Cos@xD + Sin@xDL== 2
DSolve@y ''@xD + k ^ 2 y@xD 0, y@xD, xD 88y@xD ® C@1D Cos@k xD + C@2D Sin@k xD<< DSolve@8y ''@xD + k ^ 2 y@xD 0, y@0D == 1<, y@xD, xD 88y@xD ® Cos@k xD + C@2D Sin@k xD<< DSolve@8y ''@xD + k ^ 2 y@xD 0, y@0D == 1, y '@0D 0<, y@xD, xD 88y@xD ® Cos@k xD<<
H* trigonometricke funkce *L Cos@x + yD Cos@x + yD vysl = TrigExpand@Cos@x + yDD Cos@xD Cos@yD - Sin@xD Sin@yD TrigFactor@vyslD Cos@x + yD f = Hn1 Cos@Α 2D - n2 Cos@Α 1DL Hn1 Cos@Α 2D + n2 Cos@Α 1DL - n2 Cos@Α 1D + n1 Cos@Α 2D n2 Cos@Α 1D + n1 Cos@Α 2D n1 = n2 Sin@Α 2D Sin@Α 1D n2 Csc@Α 1D Sin@Α 2D f2 = Simplify@fD - Sin@2 Α 1D + Sin@2 Α 2D Sin@2 Α 1D + Sin@2 Α 2D TrigFactor@f2D - Cot@Α 1+ Α 2D Tan@Α 1- Α 2D TrigToExp@Sin@xDD 1 2
ä ã-ä x -
1
ä ãä x
2
9
math_tisk.nb
ExpToTrig@Exp@I xDD Cos@xD + ä Sin@xD
H* rozvoj v radu *L Series@Exp@xD, 8x, 0, 3
x3
1+x+
+ 2
+ O@xD4
6
Series@Sqrt@1 + xD, 8x, 0, 3
x 1+
-
x3 +
2
8
+ O@xD4
16
Series@Hn + 2L Hn + 3L, 8n, Infinity, 3
3 + n2
n
1 4 + OA E n n3 9
-
Normal@%D 9
3
1-
+ n3
1 -
n2
n
H* grafy funkci *L plot1 = Plot@Sin@xD, 8x, - 5, 5
0.5
-4
2
-2
4
-0.5
-1.0
10
math_tisk.nb
plot2 = Plot@Sin@x + 1D, 8x, - 5, 5
0.5
-4
2
-2
4
-0.5
-1.0
Show@8plot1, plot2<, AxesLabel ® 8x, y
0.5
x -4
2
-2
4
-0.5
-1.0
f = Sin@Sqrt@x ^ 2 + y ^ 2DD Sqrt@x ^ 2 + y ^ 2D SinA
x2 + y 2 E x2 + y2
11
math_tisk.nb
Plot3D@f, 8x, - 30, 30<, 8y, - 30, 30
0.1 0.0
20
-0.1 0 -20 0 -20 20
Plot3D@f, 8x, - 30, 30<, 8y, - 30, 30<, PlotRange ® 8- 0.2, 1
1.0
0.5 20 0.0 0 -20 0 -20 20
12
math_tisk.nb
Plot3D@f, 8x, - 30, 30<, 8y, - 30, 30<, PlotRange ® 8- 0.1, 1<, PlotPoints ® 50D
1.0
0.5 20 0.0 0 -20 0 -20 20
DensityPlot@f, 8x, - 20, 20<, 8y, - 20, 20<, PlotRange ® 8- 0.22, 0.7<, PlotPoints ® 100D 20
10
0
-10
-20 -20
-10
0
10
20
H* vizualizace dat *L a = Table@Sin@xD, 8x, - 5, 5, 0.1
13
math_tisk.nb
ListPlot@aD 1.0
0.5
20
40
60
80
100
-0.5
-1.0
a = Table@8x, Sin@xD<, 8x, - 5, 5, 0.1
0.5
-4
2
-2
4
-0.5
-1.0
f SinA
x2 + y 2 E x2 + y2
a = Table@Sin@xD, 8x, 0, 5, 0.5
0.5
-0.5
-1.0
a = Table@f, 8x, - 10.001, 10<, 8y, - 10.001, 10
14
math_tisk.nb
ListPlot3D@a, InterpolationOrder ® 0D
0.4 0.2
20
0.0 15 -0.2 10
5 10
5 15 20
ListPlot3D@aD
0.4 0.2
20
0.0 15 -0.2 10
5 10
5 15 20
15
math_tisk.nb
ListDensityPlot@a, InterpolationOrder ® 0D
20
15
10
5
5
10
15
20
ListDensityPlot@a, InterpolationOrder ® 3D
20
15
10
5
5
10
15
20
a = Table@Exp@Hi - jLD, 8i, 1, 4<, 8j, 1, 4
1 ,
991, ã
1 ,
ã2
1 =, 9ã, 1,
ã3
1 ,
ã
=, 9ã2 , ã, 1,
ã2
1
=, 9ã3 , ã2 , ã, 1==
ã
16
math_tisk.nb
BarChart3D@a, ChartLayout ® "Grid", Boxed ® False, ChartStyle ® Directive@EdgeForm@BlackD, BlueD, ChartLabels ® 81, 2, 3, 4<, Method ® 8"Canvas" ® None<, Ticks ® None, ViewPoint ® 8- 10, - 30, 15<, PlotRange ® 8- Max@aD, Max@aD
4 3 2 1
H* cteni dat ze souboru *L data = ReadList@"data.txt", Number, RecordLists ® TrueD; Max@dataD 232. Min@dataD 35.
17
math_tisk.nb
obr1 = ListDensityPlot@data, PlotRange ® All, Frame ® False, InterpolationOrder ® 0D
obr2 = ListPlot3D@data, Mesh ® False, Ticks ® FalseD
H* export grafiky do souboru *L Export@"obrazek1.jpg", obr1D obrazek1.jpg
18
math_tisk.nb
Export@"obrazek2.jpg", obr2D obrazek2.jpg
H* cviceni *L H* prvocisla *L prvocisla@max_D := H Print@2D; For@i = 3, i £ max, i += 2, prvoc = True; For@del = 3, del £ Sqrt@iD, del += 2, If@Divisible@i, delD, prvoc = FalseD D If@prvoc, Print@iDD DL prvocisla@100D
19
math_tisk.nb
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97 prvocisla@max_D := H list = 82<; For@i = 3, i £ max, i += 2, prvoc = True; For@del = 3, del £ Sqrt@iD, del += 2, If@Divisible@i, delD, prvoc = FalseD D If@prvoc, list = Append@list, iDD D; listL prvocisla@100D 82, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97<
H* statistika *L d = ReadList@"statistika.txt", NumberD;
20
math_tisk.nb
pocet = Length@dD 10 000 Mean@dD N 3.9678 Variance@dD N 3.94736 plot1 = Histogram@dD
2000
1500
1000
500
5
10
15
p = Table@pocet PDF@PoissonDistribution@4D, iD, 8i, 0, 15
1500
1000
500
5
10
15
21
math_tisk.nb
Show@8plot1, plot2
2000
1500
1000
500
5
10
15
p2 = Table@8i + 0.5, pocet PDF@PoissonDistribution@4D, iD<, 8i, 0, 15
1500
1000
500
5
10
15
Show@8plot1, plot3
2000
1500
1000
500
5
10
15
H* difrakce *L z =.
22
math_tisk.nb
u = Integrate@Exp@- I Hx - ΞL ^ 2 zD, 8Ξ, - 1, 1
I
1 2
-
ä M 2
Π 2
KErfA
H-1L14 H-1+xL
E - ErfA
H-1L14 H1+xL
z
EO
z
z z = 0.01 0.01 Plot@Abs@uD, 8x, - 2, 2<, AxesOrigin ® 80, 0<, PlotRange ® AllD 20
15
10
5
-2
1
-1
2
z = 0.3 0.3 Plot@Abs@uD, 8x, - 10 z, 10 z<, AxesOrigin ® 80, 0<, PlotRange ® AllD 4
3
2
1
-3
-2
-1
1
2
3
z = 100 100
23
math_tisk.nb
Plot@Abs@uD, 8x, - 10 z, 10 z<, AxesOrigin ® 80, 0<, PlotRange ® AllD 0.020
0.015
0.010
0.005
-1000
500
-500
1000
z =. u
I
1 2
-
ä M 2
Π 2
KErfA
H-1L14 H-1+xL
E - ErfA
H-1L14 H1+xL
z
EO
z
z
I
1 2
-
ä M 2
Π 2
KErfB
H-1L14 H-1+xL
F - ErfB
z
fresnel@x_, z_D := z
24
math_tisk.nb
H-1L14 H1+xL z
FO
Animate@Plot@Abs@fresnel@x, zDD, 8x, - 10 z, 10 z<, AxesOrigin ® 80, 0<, PlotRange ® AllD, 8z, 0.1, 1<, AnimationRunning ® FalseD
z
1.0
0.5
-1.0
0.5
-0.5
1.0
-0.5
-1.0
25
math_tisk.nb
