伊達 >> 応用数学2 >> 課題一覧 >> 課題1

課題1． C言語バージョン （2010年11月25日作成）   

目的

調べること

• いくつか簡単な関数について，そのフーリエ級数展開がどうなるか， 手計算か，教科書などで調べておく．
• フーリエ級数の，項の数を増やすにしたがい， もとの関数を精度よく表現できることを確認する．
• 近似の精度のよさを測ってみる．

  近似の良さは，もとの関数と，フーリエ級数で作った関数が どのくらい類似しているか，距離を計ればわかる． ヒント： 「関数をベクトルとして見る」

例： y(t)=|t| のフーリエ級数展開 （教科書 p.70）

• C 言語でデータを生成するプログラムを書く fourier003.c

  #include <stdio.h>
  #include <stdlib.h>
  #include <math.h>
  
  int N = 20; // どこまで足すか
  int SHOW_EVERY_N_STEP = 5;
  
  double start_t = -3.0*M_PI; // 横軸はじまり
  double end_t =  3.0*M_PI; // 横軸おわり
  double step = 0.01; 
  
  int main(int argc, char *argv[])
  {
      int i;
      double pi = M_PI;
      double t, y;
      char filename[256] = "2010am2.dat";
      FILE *fp;
      
      if((fp = fopen(filename, "wt")) == NULL){
       printf("NOT open file (%s)\n", filename);
       exit(1);
      }
  
      for ( t = start_t; t < end_t; t=t+step ){
        y = pi/2.0;   // 平均値． a0/2 項
        fprintf(fp, "%.2lf \t",t);
        for ( i = 1; i <= N; i++ ){  // N 項まで足していく
           y = y - 4.0/pi/(2.0*i-1.0)/(2.0*i-1.0)*cos( (2.0*i-1.0)*t );
           if (i % SHOW_EVERY_N_STEP == 0){
              fprintf(fp, "%.5lf ",y);  // ある項数ごとに値をプリント
           }
        }
        fprintf(fp, "\n");
      }
      fclose(fp);
  }
   
  

• コンパイルして走らせる．

  % gcc fourier003.c -lm
  % a.out
  

• 2010am2.dat というファイルが新しくできているはず．その中身をみてみる．

  % less 2010am2.dat | more
  
  -9.42   3.07814 3.10979 3.12038 3.12568
  -9.41   3.07782 3.10915 3.11943 3.12441
  -9.40   3.07687 3.10725 3.11660 3.12068
  -9.39   3.07528 3.10412 3.11198 3.11467
  -9.38   3.07307 3.09978 3.10569 3.10668
  -9.37   3.07024 3.09430 3.09793 3.09712
  -9.36   3.06679 3.08776 3.08890 3.08639
  -9.35   3.06275 3.08022 3.07884 3.07494
  -9.34   3.05812 3.07179 3.06799 3.06319
  -9.33   3.05292 3.06256 3.05660 3.05147
  -9.32   3.04717 3.05265 3.04491 3.04006
  -9.31   3.04088 3.04215 3.03312 3.02912
  -9.30   3.03407 3.03119 3.02142 3.01873
  -9.29   3.02676 3.01986 3.00996 3.00886
  -9.28   3.01899 3.00828 2.99883 2.99942
  .....................
  

• 各行に何の値が書かれているかは，プログラムを読んで下さい．
• gnuplot を使って表示してよう．

  % gnuplot
  
  
          G N U P L O T
          Version 4.2 patchlevel 5
          last modified Mar 2009
          System: Linux 2.6.31-22-generic
  
  Terminal type set to 'wxt'
  
  gnuplot> plot "2010am2.dat" using 1:2 with lines 1, "2010am2.dat" using 1:3 with lines 2
   

  で画面に表示される． LaTeX でレポートを書く際には，グラフの eps ファイルを作成しておけばいい． eps ファイルの作成の仕方は

  gnuplot> set terminal postscript "Helvetica" 20 color eps enhanced
  gnuplot> set output "foo.eps"
  gnuplot> set xlabel "t"
  gnuplot> set ylabel "f(t)"
  gnuplot> plot "2010am2.dat" using 1:2 with lines 1, "2010am2.dat" using 1:3 with lines 2
  

• 毎回，これをやると大変なので，グラフ作成用の スクリプトをあらかじめ書いておくとよい．
• gp-script． LaTeX に張り付けるための eps ファイルを生成

  set terminal postscript "Helvetica" 20 color eps enhanced
  set xlabel "t"
  set ylabel "f(t)"
  plot [:][:] "2010am2.dat" using 1:2 title "n=5" with lines 1,\
  "2010am2.dat" using 1:2 title "n=10" with lines 2,\
  "2010am2.dat" using 1:3 title "n=15" with lines 3,\
  "2010am2.dat" using 1:3 title "n=20" with lines 4
  

  このファイルを用意し

  % gnuplot gp-script > foo.eps
  

  とすると，eps ファイルが作成される．

• セーブした postscript ファイルを画面で見る．

  % evince foo.eps   # グラフを見る
  

• レポートに図を張り付ける

• 過去に入力したコマンドをいちいち繰り返し入力する必要はない． コンピュータは過去の履歴を記憶しており， 上向きの矢印キーを押せば， 最近使用したコマンドから順に再利用できる形で表示される．

そのほか

• LaTeX で楽にレポートを作成する方法
• gnuplot入門
• GNU Octave 日本語マニュアル     プロット