少し修正。 最初に一文字分読み込むことでデータサイズが0かどうかに関連するちょっとした無駄を無くす。 混乱を避けるためbool_t型にTRUEはない。FALSEであるか否かでブール値を表す。 #include <stdio.h> typedef enum { FALSE = 0 } bool_t; int main(void) { unsi</stdio.h>…

行数とは一連の文字列データに含まれる改行の数とする。 ただし文字列の最後が改行でない場合は改行の数に1を加えたものを行数とする。 また一文字もない場合は行数は0であるとする。 標準入力からデータを読み込み行数を数え表示するプログラム。 #include <stdio.h></stdio.h>…

Javaからコマンドを与えgnuplotにチャートを描かせてJava側でそれを表示する。 コマンドとチャートデータの授受は標準入出力を使用する。 チャートデータはPNG形式を利用する。 ついでに、標準エラーから得られる-オプション付きgnuplotからのメッセージも表…

滑らかなものばかりだったので滑らかさのないものを扱ってみる。 両側を壁で囲まれた一次元の質点の運動について、 壁で瞬間的に衝突、反射するとすれば運動方程式からの数値解析が難しくなるので、 壁がばねで支えられているとして、その復元力をもって有限…

CやJavaに慣れるとFORTRANでの接尾辞なしの浮動小数点型リテラルが単精度であることを忘れてしまう。 前項でのLorenz方程式の初期値問題で与えた初期値 double precision, parameter :: x0(3) = (/1d-1, 1d-1, 1d-1/) ! initial coordinate at t0 が double …

有名どころでLorenz方程式を解いてみる。 方程式のパラメータはLorenzが最初の論文で示したものを使用する。 program ode use runge_kutta_4 implicit none double precision, parameter :: a = 10, b = 28, c = 8d0/3 ! Lorenz's parameters double precisi…

前項の補足。カテゴリー的にはmathsというよりphysicsだとは思うが。 \documentclass{jsarticle} \begin{document} \subsection*{万有引力下の物体の運動方程式} 質量$M$の物体との間に働く万有引力を受けて運動する質量$m$の物体について考える。 $m \ll M$…

太陽を周回する地球の運動について見てみる。 解くべき運動方程式は時刻tでの軌道面における地球の直交座標を(x(t),y(t))として、 である。太陽は原点から動かないものとする。GMは万有引力定数と太陽質量の積である。 program ode use runge_kutta_4 implic…

RK4モジュールを使って落体の運動方程式。 せっかく計算本体を別モジュールにしたにも関わらず記述量が変わっていないのがなんとも。 program ode use runge_kutta_4 implicit none double precision, parameter :: g = 9.80665 ! gravitational acceleratio…

RK4を使って一階連立微分方程式を解くためのモジュール。 サブルーチンset_funcでのを設定し、 関数nextで時刻tの値から時刻t+hの値を算出する。 FORTRAN 2003以降でないとコンパイルできないと思う。 module runge_kutta_4 implicit none interface functio…

落体の運動方程式をRK4により数値的に解くFORTRANプログラム。 Euler法によるプログラムのうちy(t)とv(t)からy(t+h)とv(t+h)を求める漸化式部分をRK4に変更しただけ。 同じ計算をnextyとnextvで行っているのでおいしくないプログラムになっている。 もっとま…

\documentclass{jsarticle} \begin{document} \subsection*{落体の運動方程式をRK4で解く} 先述の落体の運動方程式をEuler法の替わりにRK4で解いてみる。 落下距離と落下速度の連立微分方程式はRK4では、 \begin{eqnarray*} y_{n+1} &=& y_n + \frac{h}{6}\l…

落体の運動方程式をEuler法により数値的に解くFORTRANプログラム。 program ode implicit none double precision, parameter :: g = 9.80665 ! gravitational acceleration double precision, parameter :: t0 = 0 ! initial t double precision, parameter …

\documentclass{jsarticle} \begin{document} \subsection*{落体の運動方程式をEuler法で解く} 鉛直下方向を鉛直方向の位置の変数$y$の正の方向とする。 時刻$t_0$で位置$y_0$、速度$v_0$であるような自由落下する質点の運動方程式は、 \[ m\frac{\mathrm{d}…

\documentclass{jsarticle} \begin{document} \subsection*{連立一階常微分方程式をEuler法で解く} $m$を$1 \le m \le M$の整数とし、 $M$個の変数$x_m$が$t$を独立変数とする関数であるとする。 これらが以下の微分方程式を満たすとする。 \[ \frac{{\mathr…

定番のclassical 4th-order Runge-Kutta法、略してRK4を使ってみる。 ...snip double precision function nextx(t, x, h) double precision, intent(in) :: t, x, h double precision :: k1, k2, k3, k4 k1 = f(t, x) k2 = f(t + h / 2, x + h / 2 * k1) k3 …

Heun法を使ってみる。 ...snip double precision function nextx(t, x, h) double precision, intent(in) :: t, x, h double precision :: k1, k2 k1 = f(t, x) k2 = f(t + h, x + h * k1) nextx = x + h / 2 * (k1 + k2) end function nextx ...snip 計算し…

Euler法よりもマシなはずの中点法で計算してみる。 Euler法との違いは漸化式のの定義部分だけである。 ...snip double precision function nextx(t, x, h) double precision, intent(in) :: t, x, h double precision :: k1, k2 k1 = f(t, x) k2 = f(t + h /…

関数についての常微分方程式の初期値問題 の解析解はである。 これをEuler法で数値解析するプログラムをFORTRAN、古流書式でなく90/95で書いてみる。 微分方程式の初期値問題が の形の場合、刻み幅をとし、 と定義すると、Euler法は、 と書けるので、素直に…

GTK+3用に書き直してみた。 以前書いたように、 廃されたexpose-eventシグナルの代わりに、 drawシグナルに対応させなければならない。 ...snip gboolean on_expose(GtkWidget *widget, cairo_t *cr, gpointer data); ...snip gtk_widget_modify_bg(da, GTK_…

cairotest5.cでは背景の色を白にするために、 make_offscreen()内でイメージサーフェスを矩形で塗り潰してから、 draw_hello()へサーフェスを渡すという面倒なことを行っている。 しかし、これはもっと簡単に実現できる。 gtk_widget_modify_bg()でGtkDrawin…

Win32API版と同じようにオフスクリーンバッファを利用する形にした。 そして、draw_hello()をそのままの形で使用できるようにした。 また、Win32API版に合わせて背景を白色に変更した。 cairotest5.c #include <gtk/gtk.h> #define UNUSED_ARGUMENT(x) (void)(x) #defin</gtk/gtk.h>…

Win32APIの代わりにGTK+2をトップレベルウインドウの作成に使う。 フレームワークを使うとコードがずいぶん短くなる。 ただし、処理的には同じことをしているのだが、 draw_hello()を直接には使っていない。 GTK+3ではexpose-eventシグナルがdrawシグナルに…

Win32APIを直接利用してcairoで描画する領域を用意してみる。 cairoにはWin32サーフェスが用意されており、 描画対象のデバイスコンテクストのハンドルを渡せばcairoサーフェスとして機能する。 Win32サーフェスを使うときはcairo-win32.hもインクルードする…

cairoの公式ページのFAQにあるPNGファイル出力のような形に変更するのも簡単である。 --- cairotest1.c # +++ cairotest2.c # @@ -1,13 +1,13 @@ #include <cairo.h> -#include <cairo-pdf.h> void draw_hello(cairo_surface_t *surf); int main (void) { - cairo_surface_t *surf </cairo-pdf.h></cairo.h>…

Windowsでcairoを利用したソースをコンパイルする時に面倒なのは、 コンパイルオプションを生成するためのpkg-configの利用だと思う。 UNIXライクシェルが使えるなら多分サブシェルを利用したりバッククォート記法が使えると思うので問題ないが、 cmd.exeで…

描画ライブラリのcairoはGTK+等のcairoを想定したフレームワークから使うと便利だが、 そういったライブラリを介さずに使ってもそれほど変わりがあるわけではない。 コードの既述量の違いはcairoを取り巻く部分であり、 cairoに直接関わっている部分は同じよ…

の導関数はなので、 x=0のときy'=0となり、ニュートン・ラフソン法においてはx=0は初期値として不適である。 ところが2階中心差分による微分係数ではでであり、 誤って不適切なパラメータが与えられると初期値x=0でも不幸にも正しい解に収束してしまったりす…

はさみうち法でやってみる。 PROGRAM INTERS IMPLICIT NONE DOUBLE PRECISION ERR,HIGH0,LOW0,Y1,Y2,X,Y,HIGH,LOW,MID PARAMETER(ERR=1.0D-10,HIGH0=1.0,LOW0=0.0) Y1(X)=X*TAN(X) Y2(X)=SQRT(1-X**2) HIGH=HIGH0 LOW=LOW0 DO MID=(HIGH+LOW)/2 Y=Y2(MID)-Y1…

導関数を与えて微分係数を算出する方法。 PROGRAM INTERS IMPLICIT NONE DOUBLE PRECISION Y,YD,Y1,Y2,X,XN,E,X0 PARAMETER(E=1.0D-12,X0=0.9) Y1(X)=X*TAN(X) Y2(X)=SQRT(1-X**2) Y(X)=Y1(X)-Y2(X) YD(X)=(TAN(X)+X/COS(X)**2)-(-X/Y2(X)) X=X0 DO XN=X-Y(X…