クランク・ニコルソン法

 偏微分方程式を解こう第三弾としてクランクニコルソン法に挑戦した。

これは今までやってきた二つの方法とは毛色が違って計算が面倒くさい。

各時間ステップごとに空間分割数元の連立方程式を解く必要がある。

連立方程式を解く部分はGSLに頼った。

ソースコードはまたそのうちアップする。


あと、Lax法の最低次誤差項が2階微分であること、

Lax-Wendroff法は3階微分であることを確かめた。
ジャジャガッチ | コンピュータ | 22:14 | comments(0) | trackbacks(0) |

グラフをアニメ化する方法

 ぱらぱら漫画の要領でアニメ化する。

まず、各コマになるデータファイルを用意する(0.dat,1.dat・・・)。

これをgnuplotでグラフ化する。

set yrange [0:1]
set terminal gif

set output "000.gif"
plot "0.dat"
set output "001.gif"
plot "0.dat"



とすれば000.gif、001.gif・・・と連番になったgifファイルを得られる。

これをGiamというフリーソフトでアニメかしたいのだが、なぜかこのままだと読めない。

IrfanViewなどを使ってもう一度gifとして保存しなおす。

File=>Batch Conversion・・・で出てきたウィンドウにファイルを

ドラッグアンドドロップし、Start Batchボタンを押す。

このようにして得られたgifをGiamに突っ込むとぱらぱらアニメが完成する。

参考記事(ジャジャガッチの勉強ノート)
ジャジャガッチ | コンピュータ | 13:44 | comments(0) | trackbacks(0) |

関数ポインタの使いどころ

 長いこと関数ポインタは無視していたのだが、ちょっとわかってきた。

例を挙げよう。

#include <stdio.h>
#include <stdlib.h>

void f()
{
    printf("Hello.¥n");
}

void g()
{
    printf("hello.¥n");
}

int main()
{
    int i;
    void (*p[2])();
    p[0] = f;
    p[1] = g;

    for(i=0;i<2;i++) (p[i])();
   
    return 0;
}

if文を使わなくても関数の使い分けができる。

もうひとつは構造体とのコラボである。構造体は本来関数をメンバにもてないが、

関数ポインタを使うことで可能になる。

#include <stdio.h>
#include <stdlib.h>

struct TEST
{
    int a;
    void (*p)();
};

void f()
{
    printf("Hello.¥n");
}

int main()
{
    struct TEST test = {1,f};

    test.p();
   

    return 0;
}

GSLでこういう使い方をされているのを見て知った。
ジャジャガッチ | C/C++ | 22:00 | comments(0) | trackbacks(0) |

構造体の使いどころ

関数の引数として渡すのが便利だと思った。
 
#include <stdio.h>

int  f(int a,int b,int c)
{
    return a+b+c;
}

int main()
{
    f(1,2,3);

    return 0;
}

このような書き方をすると、例えば関数の仕様を変更してf(a,b,c,d)と4つの引数を

受け取るようにしたとき、関数の宣言、定義だけでなく、関数を使用しているところ

全てを変更する必要がある。f(1,2,3)=>f(1,2,3,4)みたいに。構造体を使うとすっきり

書ける。

#include <stdio.h>
#include <stdlib.h>

struct ABC
{
    int a,b,c;
};

int  f(struct ABC abc)
{
    return abc.a+abc.b+abc.c;
}

int main()
{
    struct ABC abc = {1,2,3};
   
    f(abc);

    return 0;
}

これだと仕様変更にも構造体の変更で対応できる。

ポインタを引数としても同様な形に書けるが、配列の形にまとめるのは違和感が

ある変数群は構造体の形でまとめておけばよい。

ここでは値渡ししているが参照渡しも出来る(構造体ポインタ)。
ジャジャガッチ | C/C++ | 21:46 | comments(0) | trackbacks(0) |

ポインタの使いどころ(2)

 ポインタの使いどころ

2.でっかい配列を使いたいとき

C言語ではメモリ領域が二つに分かれている。スタック領域とヒープ領域だ。

普通の変数や配列はスタック領域に置かれる。

スタック領域はあまり大きくない。うちの環境だと次のプログラムはエラーが出る。

#include <stdio.h>

int main()
{
    double a[500000];
   
    return 0;
}

スタック領域には格納しきれないのだ。次のプログラムは正しく動く。
#include <stdio.h>
#include <stdlib.h>

int main()
{
    double *a;

    a = (double*)malloc(sizeof(double)*500000);

    free(a);
    return 0;
}

mallocを使うことでヒープ領域にメモリを確保できる。ヒープは大きい。

値へのアクセスは配列のように出来る。
#include <stdio.h>
#include <stdlib.h>

int main()
{
    double *a;
    int i;

    a = (double*)malloc(sizeof(double)*5);

    for(i=0;i<5;i++) a[i] = i;

    free(a);
    return 0;
}

この応用で、可変長の配列が必要なときにもポインタが役立つ。
C言語の配列はコンパイル時に大きさが決まっていないといけない。
これは不便だ。例えばテキストファイルを読み込んで配列に格納することを
考えると配列の大きさはテキストファイルの大きさに合わせて変えたい。
このようなときmallocを使えばテキストファイルの大きさに合わせて
メモリの確保をすることが出来る。
ジャジャガッチ | C/C++ | 21:36 | comments(0) | trackbacks(0) |

ポインタの使いどころ

ポインタの使いどころ

1.参照渡しをしたいとき
 
#include <stdio.h>

void f(int a)
{
    a = 2;
}

int main()
{
    int a = 1;
    f(a);
    printf("%d",a);
    
    return 0;
}

f関数にaを渡してa=2にしているのでコンソールには2が表示されて欲しいのだが、実際は1が表示される。
C言語の関数は値渡ししか出来ないので、新たにint bが生成され、b=a;b=2;が実行されているのと同じで、main関数のaには影響しない。

main関数のaの値を変えたいときにポインタを使う。

#include <stdio.h>

void f(int *b)
{
    *b = 2;
}

int main()
{
    int a = 1;
    int *b;

    b = &a;

    f(b);
    printf("%d",a);
   
    return 0;
}
このようにするとちゃんとmain関数のaの値が2になる。アドレスを渡しているのでこういうことが出来る。
この例だと、わざわざポインタを使わなくて
#include <stdio.h>

int f(int a)
{
    a = 2;
    return a;
}

int main()
{
    int a = 1;
   
    a = f(a);
    printf("%d",a);
   
    return 0;
}
こうしてもいいんだけど、変更したい数が多くなるとこうはいかなくなる。それがポインタを使うと・・・
#include <stdio.h>

void f(int *a)
{
    int i;

    for(i=0;i<5;i++)a[i] = 2;
}

int main()
{
    int a[5];
    int i;
   
    f(a);

    for(i=0;i<5;i++) printf("%d",a[i]);
   
    return 0;
}

厳密に言えばこれは配列だけど。

C言語を学び始めたときは関数の引数は入力パラメータというイメージが強かったが、ポインタを少し理解すると、イメージが変わった。上の例のように計算結果格納用の配列を用意して空っぽのまま関数に渡すこともできるのだ。
次回に続く。
ジャジャガッチ | C/C++ | 20:54 | comments(0) | trackbacks(0) |

C言語の基礎(ポインタ)

 C言語を独学で学び始めてから数年経つが、構造体と関数ポインタの便利さが

もうひとつわからなかった。

今日あるサンプルプログラムを読んでいて少し便利さがわかったのでまとめておく。

ポインタから話を始める。

ポインタは初心者にはわかりづらく、その意味で非常に悪名高い。

実際僕も例に漏れずつまづいた。

メモリがどうのこうのというのはわかるし、文法もわかるのだが、

便利さがよくわからなかった。

当時の自分のプログラムではポインタを使う必要性もなかったし仕方なかったと思う。

必要性がないと便利さがよく理解できないものだ。

人の作ったライブラリを使おうと思うとポインタとかばしばし使ってるから

勉強になります。というかポインタわからないと使えない。

前置きが長くなったけどまずはポインタの便利さから。

次回に続く。
ジャジャガッチ | C/C++ | 20:31 | comments(0) | trackbacks(0) |

一次元移流方程式(2)

 夜ベッドで寝入る前に移流方程式について考えていたんだけれど、

密度の塊は時間が経っても分散しないんじゃないかと思った。

昨日の計算結果だと少しずつ分散していっていたけれども。

で、今日調べたらLax法はそういう性質があるらしい。

散逸誤差とか数値拡散とかいうらしい。

差分方程式を逆にテイラー展開を使って微分形に戻してやると確かにそれらしい項がある。

誤差項が偶数回微分→散逸、奇数回微分→振動(位相誤差)ということらしい。

ここらへんきちんと勉強したい。

あと、フォンノイマン安定性解析でスキームを評価したけれど、きちんと解を予測しておかないと

まずいと思った。発散しなければいいというものでもない。

誤差で減衰していくのはまずいし、発散するのが物理的に正しいこともあるんでねえか。

Lax法は散逸するということで、振動するLax-wendroff法も試してみた。

確かに振動している。そのうちソース載せる。
ジャジャガッチ | コンピュータ | 22:10 | comments(0) | trackbacks(0) |

一次元移流方程式

 円周率の計算は意外と面倒くさそうなのでやめ。

また気が向いたら勉強する。

もっと手軽に出来そうなものはないかと探してみた。

数値計算ライブラリGSLのマニュアルを読んでいたら

偏微分方程式を解く機能がないことに気が付いた。

この先偏微分方程式を数値的に解く機会がないとも言えないし、

勉強がてらやってみることにした。

とりあえず数値計算分野で有名な本「ニューメリカルレシピ」を引っ張り出してきた。

偏微分方程式の章を少し読んだ。

最初に一次元移流方程式というものを例にとっていたので、まずこれをやってみた。

ある一次元の管を考えて、一定速度で水が流れているとする。

そこに何か物質を入れたときにどのように密度変化するかを記述する方程式だ。

物質同士と水同士の相互作用は考慮されていない。

Lax法という方法で解いた。

時刻ごとに密度分布ファイルを生成し(10000ファイル!)、gnuplotで順次読み込んで

アニメにした。それっぽいアニメになったので多分合ってる。

ジャジャガッチの勉強ノートの方にそのうちアニメとか載せるかも。
続きを読む >>
ジャジャガッチ | コンピュータ | 22:24 | comments(0) | trackbacks(0) |

円周率10000桁

  exflibを使って計算してみた。マチンの公式というやつを使った。

exflibには多倍長でatanを計算する関数があるのですごく簡単。

#define PRECISION 10100
#include <iostream>
#include "exflib/exfloat.h"

int main()
{
    exfloat a = "1/5";
    exfloat b = "1/239";
    exfloat c = 4;

    std::cout << std::setprecision(10000) << c*(c*atan(a)-atan(b)) << std::endl;
   
    return 0;
}

計算結果を載せようと思ったんだけど10000桁を適当に改行するのが面倒なのでやめ。
ジャジャガッチ | コンピュータ | 20:59 | comments(0) | trackbacks(0) |
1/2PAGES | >> |

08
--
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
--
>>
<<
--
PR
RECOMMEND
RECENT COMMENT
MOBILE
qrcode
OTHERS
Since 2013/09/17
LATEST ENTRY
CATEGORY
ARCHIVE
LINKS
PROFILE
SEARCH