行列の式の級数展開

フォーラム(掲示板)ルール
フォーラム(掲示板)ルールはこちら  ※コードを貼り付ける場合は [code][/code] で囲って下さい。詳しくはこちら
いその

行列の式の級数展開

#1

投稿記事 by いその » 10年前

私は大学の課題でn×n行列Aが与えられた時,exp(At)の値を求めるプログラムを作成する問題を解いています.
行列Aとt(スカラー・パラメータ)とEPS(終了条件)はユーザー入力です。
より、精度が高い計算を行うため、例えば10項ごとに切ったりしてグループ分けをしたりさらに,先にEPSよりも小さくなる項,つまり,終了するところを調べてから逆から足し合わせる方法などがあるらしいのですが、どうプログラムを組めば良いか分からないです。
実際に組んだプログラムは以下です。このプログラムを改善して精度を上げれるらしいのですがどこいじればいいかわからないので助けてください。
私はC言語をあまり理解できていないので実際に組んでいただけると幸いです。ソースのコメントは私自身が書いたものです。

コード:

 #include <Windows.H>  // ヘッダファイル
  #include <math.h>
  #include <stdlib.h>
  #include<stdio.h>
  
  void func(int n, double t, double eps, double a[10][10], double expat[10][10]);//関数のプロトタイプ宣言
  
  int main(void){
      int n, i, j;                                                              //i,j: for文のカウンタ  n: 行列Aのサイズ
     double eps,t;                                                              //t:スカラー・パラメータ eps:EPS(終了条件)
     double a[10][10], expat[10][10];                   //配列a: 行列A  配列expat: exp(At)
        
     
     printf("nの値を入力してください。n = ");   //行列のサイズの入力
     scanf("%d", &n);
     printf("スカラー・パラメータの値を入力してください。t = ");        //パラメータの入力
     scanf("%lf", &t);
     printf("EPSの値を入力してください。EPS = ");       //EPSの入力
     scanf("%lf", &eps);
     
     for(i = 0; i < n; i++){                                    //行列Aの要素の入力
         for(j = 0; j < n; j++){
             printf("A[%d][%d] = ", i, j);
             scanf("%lf", &a[i][j]);
         }
     }
     
     func(n, t, eps, a, expat);                                 //計算を行うfuncを呼び出す。各引き数の役割は同じ
     
     for(i = 0; i < n; i++){                                    //計算結果の出力
         for(j = 0; j < n; j++){
             printf("expAt[%d][%d] = %lf¥n", i, j, expat[i][j]);
         }
     }
     return 0;
 }
     
 void func(int n, double t, double eps, double a[10][10], double expat[
10][10]){
     int i, j, l;                                                               //i,j,l: for文のカウンタ
     int k = 1;                                                                 //第k項の計算管理のカウンタ
     double norm = eps + 1;//normは第k項の絶対値である。初期値をeps+1にすることで、while文から抜けることを防いでる
         double at[10][10], ak[10][10], s[10][10];      //配列at: 行列Aにtを乗算したもの。配列ak: exp(At)の第k項。配列s: akの値を一時的保                                 // 存しておくもの。sはexp(At)の第k-1項を表している。
     
     for(i = 0; i < n; i++){    
         for(j = 0; j < n; j++){
             s[i][j] = 0.0;
             expat[i][j] = 0.0;
             at[i][j] = a[i][j] * t;
         }
         
         s[i][i] = 1.0;                                                 //配列sと配列expatを単位行列として初期化。
                 expat[i][i] = 1.0;     
     }
     


     while(norm > eps){                                                 //終了条件であるepsよりnormが小さくなるまでループ。
         norm = 0.0;

         for(i = 0; i < n; i++){                                //for文で行列の乗算を行う。
             for(j = 0; j < n; j++){                    
                 ak[i][j] = 0.0;                                //第k項の行列の各成分を初期化。
                  
                 for(l = 0; l < n; l++){
                     ak[i][j] += at[i][l] * s[l][j];    //配列akに配列atと配列sの乗算結果を格納。行列Atとexp(At)第k-1項目の行列を乗算している。
                 }
                 ak[i][j] = ak[i][j] / k;               //最後に配列akに配列akとkを除算したものを格納。kで割ることでexp(At)の級数展 開を表現。
             }
         }
         
         for(i = 0; i < n; i++){
             for(j = 0; j < n; j++){
                 expat[i][j] += ak[i][j];                               //第k項目を配列expatに加算。
                 norm += ak[i][j] * ak[i][j];                   //normに各成分の2乗和を格納。epsより小さくなるとwhile文を抜ける。
                 s[i][j] = ak[i][j];                                    //配列sに配列akを格納する。
             }
         }
         
     k++;                                                                                               //kをインクリメントして次の項の計算に移る。
     
     }
 }

“C言語何でも質問掲示板” へ戻る