微分方程式の近似 (0)

12/10 改題・カテゴリー変更
--------------------------

よく考えると「導関数の近似(0)～(5)」はつまらなかったので削除しました．

補足： A_F ( A_Fq_n = (q_n+1 + q_n) / 2 ) を用いれば ⊿_F, A_F だけで q '(t) = -a q(t) や q "(t) = - ω² q(t) の良好な近似が得られます．（因果性による ⊿_F, ⊿_B の選択は不要）
発展： 本丸はもちろん x '_k(t) = f_k(t, x₁(t), ... , x_m(t)) のような非線形方程式です．
蛇足： 失敗作でもいいから見たいという人はこちらをどうぞ．

[1] 第16章 常微分方程式の初期値問題 — 一段法
    http://na-inet.jp/nasoft/chap16.pdf  (16.6 陽的・陰的Runge-Kutta法の比較)
[2] Special lecture about Infinite higher order Symplectic Integrator
    http://www-cms.phys.s.u-tokyo.ac.jp/~naoki/CIPINTRO/SYMP/index.html

16進-10進変換

「まず覚えるＣの文法」 P44B.h: 関数の定義 //------------------------------------- #include "P44A.h" void XtoD(char *p){ int y[N], i, j, k, d=1000*1000; //y[0]*10^0 + y[1]*10^6 + ... k=1; y[k]=y[0]=0; for(j=0; j<m; j++){ for(i=0; i<k; i++){y[i]<<=8;} y[0]+=x[j]>>8; for(i=0; i<k; i++){ y[i+1]+=y[i]/d; y[i]%=d; } for(i=0; i<k; i++){y[i]<<=8;} y[0]+=x[j]&0xFF; for(i=0; i<k; i++){ y[i+1]+=y[i]/d; y[i]%=d; } if(y[k]>0){k++; y[k]=0;} } //p[0]*10^{j+6*k} + ... k--; i=y[k]; for(j=0; i>0; j++){i/=10;} for(i=0; i<j; i++){ p[j-i-1]='0'+y[k]%10; y[k]/=10; } k--; p+=j; if(j==0){*p='0'; p++;} while(k>=0){ for(i=0; i<6; i++){ p[5-i]='0'+y[k]%10; y[k]/=10; } k--; p+=6; } //y[1]*2^{-16} + y[2]*2^{-32} + ... if(m<n){ *p='.'; p++; k=1; for(j=m; j<n; k++,j++){y[k]=x[j];} for(j=0; j<(n-m)*5; j++){ y[0]=0; for(i=0; i<k; i++){y[i]*=10;} for(i=k-1; i>0; i--){ y[i-1]+=y[i]>>16; y[i]&=0xFFFF; } *p='0'+y[0]; p++; } } *p='\0'; } //------------------------------------- P44B.c: チェック用プログラム //------------------------------------- #include "P44B.h" int main(void){ int i; char s[6*N], c; printf("s="); scanf("%s", s); DtoX(s); for(i=0; i<n; i++){ if(i==m){printf(".");} printf("%5x", x[i]); } XtoD(s); printf("\ns==%s\n", s);   return 0; } //-------------------------------------

================================================= 実数の16進-10進変換 P44A.c のプログラムで作成した16進表示の実数を10進 数で表示するための関数 void XtoD(char *p); の定義 を P44B.h に示します． ・入力データは P44A.h で宣言された配列 x[N] です． ・入力データの整数部は y[0]*100 + y[1]*106 + y[2]*1012 + ... の形に変換し,最上位の y[k] から順に10進数に変換 します． ・入力データの小数部は y[1]*2-16 + y[2]*2-32 + y[3]*2-48 + ... の形に変換し，10倍しながら整数部を出力します． ・小数部の桁数は 5(n-m)（ ≒4(n-m)log16 ) です． ちなみに //------------------------------------- #include "P44B.h" int main(void){ int i; char s[]="3.14159265358979323846" "264338327950288419716939937510" "582097494459230781640628620899"; DtoX(s); for(i=0; i<n; i++){ if(i==m){printf(".");} printf("%5x", x[i]); } XtoD(s); printf("\ns==%s\n", s);   return 0; } //------------------------------------- を実行したときの出力は 3. 243f 6a88 85a3 8d3 1319 8a2e 370 7344 a409 3822 299f 31d0 82e fa98 s==3.14159265358979323846 264338327950288419716939937510 58209749445923077821 になりました．入力10進数の小数部が80桁で16進数の 小数部が56桁しかないためで，DtoX() の while(n<m+k) を「while(n<m+k*5/4)」に修正すると 3. 243f 6a88 85a3 8d3 1319 8a2e 370 7344 a409 3822 299f 31d0 82e fa98 ec4e 6c89 44e6 s==3.14159265358979323846 264338327950288419716939937510 582097494459230781640628620898 99132 に変わります．

10進-16進変換

「まず覚えるＣの文法」 P44A.h: 関数の定義 //------------------------------------- #include <stdio.h> #define N 50 int x[N], m, n; void err(void){ printf("error\n"); exit(1); } void DtoX(char *p){ int y[N], i, k, d=1000*1000; m=1; x[0]=0; while('0'<=*p && *p<='9'){ if(m>=N){err();}else{x[m]=0;} for(i=0; i<m; i++){x[i]*=10;} x[0]+=*p-'0'; p++; for(i=0; i<m; i++){ x[i+1]+=x[i]>>16; x[i]&=0xFFFF; } if(x[m]>0){m++;} } for(i=0; i<m/2; i++){ k=x[i]; x[i]=x[m-i-1]; x[m-i-1]=k; } n=m; if(*p=='.'){ p++; i=k=0; while('0'<=*p && *p<='9'){ if(i%6==0){k++; y[k]=0;} y[k]=y[k]*10+*p-'0'; i++; p++; } while(i%6>0){y[k]*=10; i++;} while(n<m+k){ y[0]=0; for(i=k; i>0; i--){y[i]<<=8;} for(i=k; i>0; i--){ y[i-1]+=y[i]/d; y[i]=y[i]%d; } y[0]<<=8; for(i=k; i>0; i--){y[i]<<=8;} for(i=k; i>0; i--){ y[i-1]+=y[i]/d; y[i]=y[i]%d; } x[n]=y[0]; n++; if(n>=N){err();} } } } //------------------------------------- P44A.c: チェック用プログラム //------------------------------------- #include "P44A.h" int main(void){ int i; char s[6*N]; printf("s="); scanf("%s", s); DtoX(s); for(i=0; i<n; i++){ if(i==m){printf(".");} printf("%5x", x[i]); }    return 0; } //-------------------------------------

================================================= 実数の10進-16進変換（11/24本文修正） 固定小数点の10進数で表示された正の実数を16進数に 変換して表示する関数 void DtoX(char *p); を次の設計仕様で考えてみましょう． (1) 10進数は文字列で与える． (2) 16進数は整数の配列 x[N] に上位の桁から順に4桁 ずつ格納する． (3) 16進数全体の桁数を変数 n に，整数部の桁数を変 数 m に格納する． (4) 有効数字 100 桁の実数でも変換できること． 語長が無限大の仮想的な計算機での10進-16進変換を for(w=0; '0'<=*p && *p<='9'; p++){ w=w*10+*p-'0'; } if(*p=='.'){p++; fract();} としましょう．ここで fract() は小数部の処理を行う 関数です． 整数部は w を4ビットずつ区切っていけば よいのですが，逆変換も考慮した有限語長による処理 の都合上，設計仕様(2)では x[i] の下位の 16 ビット のみにデータを格納するように規定しています． fract() では，とりあえず入力を6桁ごとに区切り y[1]×10-6 + y[2]×10-12 + y[3]×10-18 + ... に対応する y[k]×10-6k を記憶しておけば,小数部の最 初の16進4桁は小数部を216倍した実数の整数部で得られ ます． 216倍した実数の小数部に対しても同様の処理を 行えば次の4桁が求められます (以下同様)． 左記のプログラムでは ・整数部は次の形で記憶しています． x[m-1]×20 + x[m-2]×216 + x[m-3]×232 + ... ・小数部の処理では10進数6桁で20ビット必要なので 216倍（y[0]に格納）を2回に分けて行っています． ・10進数の小数部の桁数が6の倍数になるように末尾 に0を詰めています． P44A1.c のプログラムで "s=" に "0.3"（"0.30" でも 同様）と入力すると 0. 4ccc としか出力されません．"0.300000000000000000000"と 入力すると出力は 0. 4ccc cccc cccc cccc になります．"4096.999999999999"のようなチェックし やすい入力で確認してください． main()の「int i; char s[6*N]; printf("s="); scanf ("%s", s);」を int i; char s[]="3.14159265358979323846" "264338327950288419716939937510" "582097494459230781640628620899"; で置換すると出力は 3. 243f 6a88 85a3 8d3 1319 8a2e 370 7344 a409 3822 299f 31d0 82e fa98 となります．（改行は便宜的に挿入） [1] 任意精度演算 - Wikipedia

無償の開発環境

記事一覧	snapshot.html

「Ｃの文法」「Ｃから見たＣ++」では描画について説明しませんでした．これは Visual Basic 等とは異なり，言語仕様にグラフィックス関連の項目がないためです．GUI コンポーネントを用いた Windows アプリケーションを作るには言語仕様とは別に開発環境とそこで使えるコンポーネントについて学ぶ必要があります．無償の開発環境としては

[1] GNUコンパイラコレクション - Wikipedia
    http://ja.wikipedia.org/wiki/GNU%E3%82%B3%E3%83%B3%E3%83%91%E3%82%A4%E3%83%A9%E3%82%B3%E3%83%AC%E3%82%AF%E3%82%B7%E3%83%A7%E3%83%B3
[2] 窓の杜 - 【NEWS】マイクロソフト、無償の開発環境“Visual Studio 2010 ...
    http://www.forest.impress.co.jp/docs/news/20100428_364484.html
[3] WideStudio/MWT Home page
    http://www.widestudio.org/ja/

等があります．上図は天才プログラマの一人である平林俊一氏が開発した WideStudio の Web ページからの引用です．

蛇足： CUI ですが使いやすかった BCC 5.5＋BCC Developer は「Borland C++ Compiler 5.5 無償ダウンロードサービス」のリンク先が無効になってしまいました（Borland Japan --> Microfocus）．

Ｃから見たＣ++ (6)

記事一覧（Ｃの文法(1)～(17),他） まず覚えるＣの文法（progC.pdf） P62A.cpp: 例外処理 //------------------------------------- #include <stdio.h> #include <exception.h> int main(void){ FILE *fp; int x[4], i; try{ fp=fopen("data.txt", "r"); if(fp==NULL){throw "backup.txt";} }catch(char *s){ fp=fopen(s, "r"); if(fp==NULL){unexpected();} printf("read from backup.txt\n"); } for(i=0; i<4; i++){ fscanf(fp, "%d", &amp;x[i]); printf("x[%d]==%dn", i, x[i]); } <span style="color: #ff0000;">  return 0; } //------------------------------------- P1B.h: Ｃの構造体 //------------------------------------- typedef struct{ int id; char name[10]; char mail[30]; } member; //-------------------------------------

Ｃ++では正常な処理ができないときの処理を例外処理 として扱うことができ，例外処理が複雑なときは try 文を使うと処理の流れが見やすくなります[4-5,6]． ・正常な処理ができない可能性がある処理を try 節に 書き，処理ができなくなったら throw 文で後続する 処理を破棄して catch 節に移ります． ・throw 文に書かれた式の値が catch 節に渡されます． ・unexpected() で強制終了できます． Ｃ++ のテンプレートについてはプログラム例は省略し て，要点だけを紹介します． P1B.h の構造体 member にコンストラクタとメンバ関数を追加してクラスに拡 張することは容易ですが，mail の代わりに別のデータ を使いたいときは別のクラスを定義しなければなりま せん．これに対して template <class T> class member{ private: int id; char name[10]; T data; public: member(int, char*, T*); //略 }; と定めれば，member<char*>，member<struct xxx> 等 を用いて複数の型に対応できます． [4-5] 例外処理 [4-6] C++例外処理プログラミング [4-7] テンプレート関数 [4-8] テンプレートクラス