introduction to modern numerical analysis

数物セミナー冬の大談話会2015 in 大阪での発表資料

Transcript

2. • • •
3. • • • • • •
4. •
5. • • •
6. • • • •
7. ! #include <iostream> int fact(int n) { if(n == 0) { return 1; } return n*fact(n-1); } int main() { std::cout<<fact(10); } fact(10) -> 10 * fact(9) -> 10 * 9 * fact(8) -> 10 * ... * 1 * 1 -> 10! = 3628800
8. ! #include <iostream> int fact(int n) { if(n == 0) { return 1; } return n*fact(n-1); } int main() { std::cout<<fact(10); } fact 0 = 1 fact n = n*fact(n-1) main = do print $ fact 10
9. ! int fact(int n) { if(n == 0) { return 1; } else if(n < 0) { return -1; } return n*fact(n-1); } fact 0 = 1 fact n | n > 0 = n*fact(n-1)
10. •
11. • • • • • •
12. • • • • •
13. • • • •
14. • •
15. • •
16. • • • • •
17. • • • • • •
18. • • • • •
19. • • •
20. , , ∈ ℤ; 1 ≤ , , ≤ 1000 , ,
21. ( ∈ ℤ) ∈ ℤ≥0 ; 1 ≤ ≤ = + +1 + ⋯ + +−1 1 ≤ ≤ − + 1
22. 0 = 0 + 1 + ⋯ + −1 1 = 1 + 2 + ⋯ + ⋮ = + + ⋯ + + ⋮ −+1 = −+1 + −+2 + ⋯ + • × − + 1
23. 0 = 0 + 1 + ⋯ + −1 1 = 1 + 2 + ⋯ + = 0 − 0 + ⋮ = −1 − −1 + ⋮ = −1 − −1 + • 2 − •
25. std::cin std::cout
26. • • •
27. • •
28. • • constexpr • #define a b • const constexpr
31. op binary decimal 1 0000 0001 1 1 << 6 0100 0000 64 78 0100 1110 78 78 << 5 1100 0000 192
32. #include <iostream> #define N 1 << 8 // 1<<8 is 256. int main() { std::cout << N * 2; } 1. 512 (256 * 2) 2. 65536 (1<<16) 3.
33. #include <iostream> #define N 1 << 8 // 1<<8 is 256. int main() { std::cout << N * 2; } // -> 116
34. #include <iostream> #define N 1 << 8 int main() { std::cout << N*2; } // -> 116 #include <iostream> int main() { std::cout<<1<<8*2; } // -> 116
35. #include <iostream> #define N 1 << 8 int main() { std::cout << N*2; } // -> 116 #include <iostream> int main() { std::cout<<1<<8*2; } // -> 116 #define
36. #include <iostream> int main() { constexpr int N = 1 << 8; std::cout << N * 2; } // -> 512
37. #define N (1 << 8)
38. •
39. auto x = 1; // x : int auto x = 1.3; // x : double auto x; // error!
40. #include <iostream> #include <vector> int main() { std::vector<int> vec = {0, 1, 2}; for(auto const& v : vec) { std::cout << v << " "; } } // 0 1 2
41. •
42. #include <iostream> int main() { std::cout << [](int x, double y) { return x + y; }(3, 2.2); }
43. #include <iostream> int main() { auto add = [](int x, double y) {return x + y;}; std::cout << add(3, 2.2); }
44. #include <iostream> int main() { auto add = [](auto x, auto y) { return x + y; }; std::cout << add(3, 2.2); }
45. • func(int x) • func(int &x) • func(const int &x)
46. #include <iostream> void func(int x) { x++; } int main() { int a = 42; func(a); std::cout << a << std::endl; // =>42 }
47. #include <iostream> void func(int &x) { x++; } int main() { int a = 42; func(a); std::cout << a << std::endl; // =>43 }
48. #include <iostream> void func(const int &x) { x++; // error! read-only! } int main() { int a = 42; func(a); std::cout << a << std::endl; }
49. std::vector<int> v(10e7); // void func(std::vector<int> x) {} // void func(const std::vector<int> &x) {}
50. • • • • • •
51. • • • • • • • • •
52. • •
53. • + −1 −1 −1 + ⋯ + 0 = • + −1 −1 −1 + ⋯ + 0 = 0 ,
54. • • • • Δ
55. ≔ lim Δ→0 + Δ − Δ Δ , ≈ + Δ − Δ + Δ = + , Δ Δ = 0 + 0 , Δ = 0 2Δ = Δ + Δ , Δ Δ, ⋯
56. • • • • Δ +1 •
57. +1 = + Δ 6 ,1 + 2,2 + 2,3 + ,4 , ∈ ℕ 0 = 0 ,1 = , ,2 = + 1 2 Δ, + 1 2 Δ,1 ,3 = + 1 2 Δ, + 1 2 Δ,2 ,4 = + Δ, + Δ,3
58. • • • • • •
59. • • • • •
60. • • • •
61. • • • • • #include <boost/numeric/odeint.hpp>
62. void ode(const T &x, T &dxdt, const double t) { dxdt = ~~~~; } void output(const T &x, const double t) { std::cout << ~~~~; } int main() { T x; boost::numeric::odeint::integrate( ode, x, 0.0, 10.0, 0.1, output); } ∶ 0 → 10, Δ = 0.1
63. using namespace boost::numeric::odeint; boost::numeric::odeint::integrate integrate
64. using namespace std; std::cout cout
65. = − = − − = − , ,
66. = − = − − = − dxdt[0] = sigma * ( x[1] - x[0] ); dxdt[1] = x[0] * (rho - x[2]) - x[1]; dxdt[2] = x[0] * x[1] - beta * x[2];
67. typedef array<double, 3> T; constexpr double sigma = 10.0; constexpr double rho = 28.0; constexpr double beta = 8.0 / 3.0;
68. void ode(const T &x, T &dxdt, const double t) { dxdt[0] = sigma * (x[1] - x[0]); dxdt[1] = x[0] * (rho - x[2]) - x[1]; dxdt[2] = x[0] * x[1] - beta * x[2]; } void output(const T &x, const double t) { cout << t << '¥t' << x[0] << '¥t' << x[1] << '¥t' << x[2] << endl; } int main() { T x0 = {10.0, 1.0, 1.0}; // initial condition integrate(ode, x0, 0.0, 25.0, 0.1, output); }
70. std::ofstream ofs("lorenz.dat", std::ios::binary);
71. • • http://qiita.com/ignis_fatuus/items/74c b32815753d891a716
72. • http://qiita.com/ignis_fatuus/items/74c b32815753d891a716
73. = 1 + −1 + 1 + −1 2 , ∈ ℝ
74. typedef complex<double> T; constexpr complex<double> i(0.0, 1.0); void ode(const T &x, T &dxdt, const double t) { dxdt=(1.0+eta*i)*x+(1.0+alpha*i)*norm(x)*x; } int main() { T x(5.3, -2.0); }
75. + −1 −1 −1 + ⋯ + 0 = 0 • 1 1 = , 2 = , ⋯ , = −1 −1 1 = 2 , 2 = 3 , ⋯ , −1 = = − −1 + −2 −1 + ⋯ + 0 1
76. 2 2 + + = N ⋅ s2 m , N ⋅ s m , N m
77. 2 2 + + = = = − − +
78. • • #include <boost/units/systems/si.hpp> #include <boost/units/io.hpp> std::ostream
79. = 2 m, s m2, m/s 5.5 [m]
80. • m + m = + m • m × m = × m2 • m + kg → • Ω = Ω ⋅ m × m m2 • •
81. • [s] • [m] • [kg] • [A] • [K] • [mol] • [cd]
82. #include <boost/units/systems/si/prefixes.hpp> 1012 si::tera 102 si::hecto 10−3 si::milli 10−15 si::femto
83. using namespace boost::units; boost::units::si::length si::length
84. quantity<si::length, int> len; len = 1 * si::meters; // default : double quantity<si::length> len2(5.0*si::meter); quantity<si::mass> mass; mass = 5.0 * si::kilogram;
85. quantity<si::length> len(10*si::meter); quantity<si::length> len2(5.0*si::meter); cout << len + len2; // 15 m
86. quantity<si::length> len(10*si::meter); quantity<si::length> len2; len2 = 5.0 * si::kilo * si::meters; cout << len + len2; // 5010 m
87. quantity<si::length> len(10*si::meter); quantity<si::time> time(10*si::second); cout << len / time; // 10 m s^-1
88. quantity<si::length> len(10*si::meter); cout << len; // 10 m s^-1 cout << len.value(); // 10 cout << len; // 10 m s^-1
89. quantity<si::mass> len(5.0*si::meter);
90. quantity<si::length> len(10*si::meter); quantity<si::time> time(10*si::second); cout << len + time; // error!
91. #include <boost/units/pow.hpp> quantity<si::acceleration> acc(10*si::meter / pow<2>(si::second)); cout << acc; // 10 m s^-2
92. #include <boost/units/pow.hpp> quantity<si::acceleration> acc(10*si::meter / si::second); // error! length time2
93. • • • •
94. • • • • • http://kaworu.jpn.org/cpp/
95. • • http://qiita.com/hmito/items/483445ac0d4 2fb4428a5 • • http://qiita.com/t_uda/items/7712671389e 016d24df6 • • http://qiita.com/ignis_fatuus/items/74cb 32815753d891a716
96. • č • • http://www.slideshare.net/konn/haskell- 6-32258528 • • • •
97. • • • • •
98. • •
99. • • • http://nalab.mind.meiji.ac.jp/~mk/labo/t ext/
100. •