II 2019 2Q
A 6/11 6/18 6/25 7/2 7/9 7/16 7/23 B 6/12 6/19 6/24 7/3 7/10 7/17 7/24
x = 0 dv(t) dt = g Z t2 t 1 dv(t) dt dt = Z t2 t 1 gdt g v(t 2 ) = v(t 1 ) + g(t 2 t 1 ) v v(t) x g(t 2 t 1 ) t 1 t 2 t
v(t 2 ) = v(t 1 ) + g(t 2 t 1 ) = + x ( ) v V T g dv dt = g v(t) = mg 1 m v(t) e m t ϒ: t 1 t 2 t
(differential equation) h i v(t + t) = v(t) + g m v(t) v t (time step size) V T v(t) g t 1 t t t t 2 t
h i v(t + t) = v(t) + g m v(t) t mg v(t + t) = mg v(t) + apple 1 mg v(t) m t ṽ = v mg f t = m t t = t m ṽ( t + f t) = ṽ( t) + 1 ṽ( t) f t
(exact solution) v(t) = 1 t = 1 e t v 1 e t = 0.63 = 1 (time constant) t t = 0.1 1 e t = 0.63 t = 1 t
Xcode (1)
Xcode (2) mac OS > Command Line Tool Xcode Ver.
Xcode (3) Language: C++ Product name: 01FallingParticle Location:
Xcode (4) Xcode Preferences Locations Derived Data Relative
(2)Run main.cpp Build&Run (1) Build Succeeded OK (3) Console Hello OK
main.cpp #include <iostream> int main(int argc, const char * argv[]) { // insert code here... std::cout << "Hello, World!\n"; return 0; } #include <iostream> int main(int argc, const char * argv[]) { double v0 = 0.0; double dt = 0.1; double v = v0; return 0; } v 0 = 0 t = 0.1 v(0) = v 0
#include <iostream> main.cpp: int main(int argc, const char * argv[]) { double v0 = 0.0; double dt = 0.1; double v = v0; std::cout << "v: " << v << std::endl; return 0; } Build, Run Console v: 0.0 OK C printf C++ std::cout std::cout << std::endl (end of line) std standard
#include <iostream> int main(int argc, const char * argv[]) { double v0 = 0.0; double dt = 0.1; double v = v0; std::cout << "v: " << v << std::endl; v = v + (1.0 - v)*dt; std::cout << "v: " << v << std::endl; return 0; } Build, Run v
#include <iostream> int main(int argc, const char * argv[]) { double v0 = 0.0; double dt = 0.1; double v = v0; int stepend = 100; int step = 0; std::cout << "v: " << v << std::endl; while(step < stepend){ v += (1.0 - v)*dt; v +=! v = v + step++; step = step + 1 std::cout << "v(" << dt*(double)step << ") = " << v << std::endl; } return 0; } t = t step Build, Run v
gnuplot #include <iostream> int main(int argc, const char * argv[]) { double v0 = 0.0; double dt = 0.1; double v = v0; int stepend = 100; int step = 0; std::cout << dt*(double)step << " " << v << std::endl; while(step < stepend){ v += (1.0 - v)*dt; step++; std::cout << dt*(double)step << " " << v << std::endl; } return 0; }
gnuplot Terminal Terminal cd (change directory) (1)cd (2) Debug drag & drop (3)Enter pwd ( /Debug/ )
gnuplot Terminal Terminal $./01FallingParticle Terminal $./01FallingParticle > data.txt
rm: Unix Terminal pwd (print working directory) cd: (change directory) ls:./../ cd../ cp: (cp data.txt data_copy.txt ) mv: (mv data.txt../)
gnuplot Terminal gnuplot $gnuplot gnuplot >set xlabel t >set ylabel v >plot data.txt title v
gnuplot >plot data.txt title v, 1-exp(-x)
gnuplot File > Save as PDF
Xcode #include <iostream> #include <cmath> (exp ) double Output(double t, double v, bool display){ double vexact = 1.0 - exp(-t); double error = (v - vexact)/vexact; if(vexact == 0.0) error = 0.0; if(display == true) std::cout << t << ": " << v << ", " << vexact << ", " << error << std::endl; return fabs(error); }; e(t) = v(t) v ex (t) vex (t) int main(int argc, const char * argv[]) { double v0 = 0.0; double dt = 0.1; double v = v0; int stepend = 10; int step = 0; Output(dt*(double)step, v, true); while(step < stepend){ v += (1.0 - v)*dt; step++; Output(dt*(double)step, v, true); } return 0; } t = 1 Output
gnuplot gnuplot pipe gnuplot gnuplot ( ) ( )
GnuplotInterface GnuplotInterface.h &.cpp pipe www2.kobe-u.ac.jp/~hayashi /01FallingParticle/01FallingParticle/ main.cpp drag & drop
in main.cpp header #include <iostream> #include <cmath> #include GnuplotInterface.h" #include (<> ) double Output(double t, double v, bool display){...} int main(int argc, const char * argv[]) { GnuplotInterface* gnuplot = new GnuplotInterface(); double v0 = 0.0; double dt = 0.1; double v = v0; int stepend = 10; int step = 0; Output(dt*(double)step, v, true); while(step < stepend){ v += (1.0 - v)*dt; step++; Output(dt*(double)step, v, true); } delete gnuplot; return 0; } (gnuplot )
#include <iostream> #include <cmath> #include "GnuplotInterface.h" double Output(double t, double v, bool display){...}; int main(int argc, const char * argv[]) { GnuplotInterface* gnuplot = new GnuplotInterface(); gnuplot->setaxislabel("x", "t"); gnuplot->setaxislabel("y", "v"); double v0 = 0.0; double dt = 0.1; double v = v0; int stepend = 10; int step = 0; Output(dt*(double)step, v, true); } gnuplot->injection("plot '-' title 'predicted', 1-exp(-x) \n"); gnuplot->injection(dt*(double)step, v); while(step < stepend){ v += (1.0 - v)*dt; step++; Output(dt*(double)step, v, true); gnuplot->injection(dt*(double)step, v); } gnuplot->injection("e\n"); delete gnuplot; return 0; - e shift + 7 PDF copy&paste
(1) (1) I (2) t dt stepend dt = 0.1 dt = 0.01 dt = 0.001 t = 1 dt dt t=1 gnuplot Terminal x1, y1 x2, y2 gnuplot (plot data.txt ) x3, y3 set logscale x (y ) (2) 01FallingParticleImplicit v(t + t) = v(t) + (1 v(t)) t v(t + t) v(t + t) dt = 0.1
(1) I (1) t 1 1/10 1/10
(2) (2) 01FallingParticleImplicit v(t + t) = v(t) + t 1 + t #include <iostream> #include <cmath> #include "GnuplotInterface.h" void Output(double t, double v, bool display){...}; int main(int argc, const char * argv[]) {... while(step < stepend){ v = (v + dt)/(1.0 + dt); step++; Output(dt*(double)step, v, true); gnuplot->injection(dt*(double)step, v); }... }
(2) (2) v(t + t) = v(t) + (1 v(t)) t v(t + t) = v(t) + t 1 + t
t 1 v(t + t) = v(t) + 1 1! dv dt t t + 1 2! d 2 v dt 2 t t 2 + O( t 3 ) v(t + t) = v(t) + (1 v(t)) t v(t) + 1 1! dv dt t t + 1 2! d 2 v dt 2 t t 2 + O( t 3 ) = v(t) + (1 v(t)) t apple dv 1 d 2 v = 1 v(t) + dt 2! dt 2 t t t 1 + O( t 2 ) (1st-order accuracy) t 1
(3) II v(t + t) = v(t) + t 1 + t
t = n t ( n) v n = v(t) v n+1 = v(t + t) n
df dt = g(f) f n+1 t f n = g(f n k, f n ) (k > 0) (explicit scheme) (implicit scheme) f n+1 t f n = g(f n±k, f n ) (k > 0)
v n+1 = v n + (1 v n ) t v n+1 = v n + (1 v n+1 ) t t 1 Crank-Nicolson f n+1 t f n = g(f n+1 ) + g(f n ) 2
Crank-Nicolson v n+1 t v n = 1 2 (1 v n+1 ) + (1 v n ) v n+1 = vn + 1 1 2 vn t 1 + t 2 #include <iostream> #include <cmath> #include "GnuplotInterface.h" void Output(double t, double v, bool display){...}; int main(int argc, const char * argv[]) {... while(step < stepend){ v = (v + (1.0 - v*0.5)*dt)/(1.0 + dt*0.5); step++; Output(dt*(double)step, v, true); gnuplot->injection(dt*(double)step, v); }... }
Crank-Nicolson Crank-Nicolson
Adams-Bashforth Runge-Kutta Crank-Nicolson Adams-Bashforth f n+1 t f n = 3g(f n ) g(f n 1 ) 2 Runge-Kutta Runge-Kutta f (1) f n t/2 = g(f n ) f n+1 f n t = g(f (1) )
Adams-Bashforth Adams-Bashforth Runge-Kutta Runge-Kutta
t 1 t 2
(Numerical stability) t t = 0.1
gnuplot (1/4) ( ) explicit.txt implicit.txt Terminal gnuplot Terminal plot explicit.txt, implicit, 1-exp(-x) Terminal >
gnuplot (2/4) std::ofstream C fprintf (Terminal > Xcode ) ( ) ( ) int main(int argc, const char * argv[]) { std::ofstream output("data.txt");... output << dt*(double)step << " << v << " << std::endl; while(step < stepend){... step++; output << dt*(double)step << " << v << " << std::endl; }... std::cout ofstream return 0; }
gnuplot (3/4) gnuplot Xcode ( ) gnuplot plot 'data.txt' using 1:2 title 'explicit', '' using 1:3 title 'implicit' 1,2 ( data.txt ) v Crank-Nicolson C std::vector< > double v output std::vector<double> v {v0, v0}; output << dt*(double)step << " "; for(int i=0; i<v.size(); i++) output << v[i] << " "; output << std::endl;
gnuplot (3/4)... double Output(double t, double v, bool display){...}; int main(int argc, const char * argv[]) { std::ofstream output("data.txt"); ofstream GnuplotInterface* gnuplot = new GnuplotInterface(); gnuplot->setaxislabel("x", t"); gnuplot->setaxislabel("y", "v"); double v0 = 0.0; double dt = 0.1; std::vector<double> v {v0, v0}; std::vector int stepend = 10; int step = 0; output << dt*(double)step << " "; for(int i=0; i<v.size(); i++) output << v[i] << " "; output << std::endl; while(step < stepend){ for(int i=0; i<v.size(); i++){ switch (i) { case 0: v[i] += (1.0 - v[i])*dt; break; case 1: v[i] = (v[i] + dt)/(1.0 + dt); break; default: break; } } step++; output << dt*(double)step << " "; for(int i=0; i<v.size(); i++) output << v[i] << " "; output << std::endl; } gnuplot->injection("set key left top\n"); gnuplot->injection("plot 'data.txt' using 1:2 title 'explicit', '' using 1:3 title 'implicit', 1-exp(-x)\n"); for(int i=0; i<v.size(); i++) std::cout << Output(dt*(double)step, v[i], false) << ", "; delete gnuplot; gnuplot return 0; } switch
x m d2 x dt 2 = kx F = kx (k > 0) x x
dx dt = v dv dt =!2 x! = r k m x, v
dx dt = v x n+1 t x n = v n x n+1 = x n + v n t dv dt =!2 x v n+1 t v n =! 2 x n v n+1 = v n! 2 x n t x v
x(0) = x 0 = 1 v(0) = v 0 = 0! = 1 x(t) = cos(!t) T/4 v(t) =! sin(!t) T T = 2 /! T/4 10 t = T 4 1 10
#include <iostream> #include <cmath> #include "GnuplotInterface.h" int main(int argc, const char * argv[]) { GnuplotInterface* gnuplot = new GnuplotInterface(); gnuplot->setaxislabel("x", "t"); gnuplot->setaxislabel("y", "x"); } double x0 = 1.0; double v0 = 0.0; double omega = 1.0; double omega2= omega*omega; double x = x0; double v = v0; double T = 2.0*M_PI/omega; double dt = (T/4.0)*0.1; x 0 v 0!! 2 T = 2 /! gnuplot->setgraphrange("x", 0.0, T); int stepend = (int)(t/dt) + 1; int step = 0; gnuplot->injection("set key left bottom\n"); gnuplot->injection("plot '-' title 'predicted x', cos(x), 0 \n"); gnuplot->injection(dt*(double)step, x); while(step < stepend){ double xtmp = x + v*dt; double vtmp = v - omega2*x*dt; x = xtmp; v = vtmp; step++; gnuplot->injection(dt*(double)step, x); } gnuplot->injection("e\n"); delete gnuplot; return 0; 02OscillatingBody t = T 4 1 10 x n+1 = x n + v n v n+1 = v n! 2 x n t t
(x) dt 1/10
(v)
Runge-Kutta Step 1 Step 2 dx dt = v x (1) t/2 x n = v n x n+1 t x n = v (1) dv dt =!2 x v (1) t/2 v n =! 2 x n v n+1 v n t =! 2 x (1) Runge-Kutta
Runge-Kutta
Runge-Kutta f (1) = f n + f (2) = f n + t 2 gn t 2 g(1) f (3) = f n + tg (2) f n+1 = f n + t 6 (gn + 2g (1) + 2g (2) + g (3) ) 20
dx dt = v x n+1 t x n = v n+1 dv dt =!2 x v n+1 t v n =! 2 x n+1
(x)