学校の課題について

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

学校の課題について

#1

投稿記事 by sitshtr » 3年前

C言語でガウス・ザイデル法で偏微分方程式を解く課題が出されました。

ガウス・ザイデル法はなんとなく理解できたのですが、それをプログラムに持ってくる方針が立ちません、、

どのようにプログラムを組み立てていけば良いでしょうか?


よろしくお願いします!

アバター
みけCAT
記事: 6734
登録日時: 13年前
住所: 千葉県
連絡を取る:

Re: 学校の課題について

#2

投稿記事 by みけCAT » 3年前

ガウス=ザイデル法 - Wikipedia
行列の操作を多く行うようなので、行列を表す構造体を用意し、
その構造体について行列の基本演算(足し算・引き算・掛け算)を行う関数をそれぞれ用意し、
それらの関数を用いてアルゴリズムを実装するといいでしょう。
ベクトルも、n行1列の行列として扱えます。

また、まずは作るプログラムの仕様をはっきりさせるのが重要です。
まずはどのような入力(読み込む値および書式)をどこから(標準入力・コマンドライン引数・ファイル・ネットワークなど)取り、
どのような出力(書き出す値および書式)をどこへ(標準出力・ファイル・ネットワークなど)出すのかを決めましょう。
複雑な問題?マシンの性能を上げてOpenMPで殴ればいい!(死亡フラグ)

Meta3

Re: 学校の課題について

#3

投稿記事 by Meta3 » 3年前

コード:

/* 
************************************************************************
*  laplace.c:  Solution of Laplace equation with finite differences    *
*    			                                               *
*  comment: Output data is saved in 3D grid format used by gnuplot     *
************************************************************************
*/
#include <stdio.h>
#include <math.h>

#define maxx 40                         /* number of grid points */
#define maxy 40                         /* number of grid points */

main()
{
   double p[maxx][maxy], pnew, dp, dpmax, omega;
   int i, j, iter;
   
   FILE *output;			/* save data in laplace.dat */
   FILE *file;			/* save data in laplace_dp.dat */

   output = fopen("laplace.dat","w");
   file = fopen("laplace_dp.dat","w"); 

   for(i=0; i<maxx; i++)                 /* clear the array  */
   {   
      for (j=0; j<maxy; j++) p[i][j] = 0;
   }

   for(i=0; i<maxx; i++) p[i][0] = 100.0;        /* p[i][0] = 100 V */		

   for(iter=0; iter<1000; iter++)               /* iterations */
   {
      dpmax=0.0;
      for(i=1; i<(maxx-1); i++)                  /* x-direction */
      {
         for(j=1; j<(maxy-1); j++)               /* y-direction */
         {
            pnew = 0.25*(p[i+1][j]+p[i-1][j]+p[i][j+1]+p[i][j-1]);
	    dp=pnew-p[i][j];
	    if(fabs(dp) > dpmax) dpmax=fabs(dp);
            p[i][j] = pnew; 
         }
      }
      fprintf(file, "%d\t%f\n",iter,dpmax); 
   }
   
   for (i=0; i<maxx ; i++)         /* write data gnuplot 3D format */
   {
      for (j=0; j<maxy; j++) 
      {
         fprintf(output, "%f\n",p[i][j]);
      }
      fprintf(output, "\n");	  /* empty line for gnuplot */
   }
   printf("data stored in laplace.dat\n");
   fclose(output);
   fclose(file); 
}
 

Meta3

Re: 学校の課題について

#4

投稿記事 by Meta3 » 3年前

Windows10 VisualStudio2019  clang

コード:

c:\b>clang-cl c1.c

c:\b>c1.exe
data stored in laplace.dat

c:\b>
 
gnuplot でグラフ化できる

Meta3

Re: 学校の課題について

#5

投稿記事 by Meta3 » 3年前

#3
頭に  #define _CRT_SECURE_NO_WARNINGS  が必要です

Meta3

Re: 学校の課題について

#6

投稿記事 by Meta3 » 3年前

http://www.yamamo10.jp/yamamoto/lecture ... node2.html

コード:

 #define _CRT_SECURE_NO_WARNINGS

 #include <stdio.h>
 #include <math.h>
 #define N (3)               // 連立方程式の大きさ
 #define EPS (1e-15)         // 計算誤差の許容値
 
 int main(void){
   double a[N+1][N+1], x[N+1], b[N+1];
   double dx, absx, sum, new;
   int i,j;
 
   a[1][1]=3.0;  a[1][2]=2.0;  a[1][3]=1.0;  // 係数行列
   a[2][1]=1.0;  a[2][2]=4.0;  a[2][3]=1.0;
   a[3][1]=2.0;  a[3][2]=2.0;  a[3][3]=5.0;
 
   b[1]=10.0;                         // 同次項
   b[2]=12.0;
   b[3]=21.0;
 
   x[1]=0.0;                          // 近似解の初期値
   x[2]=0.0;
   x[3]=0.0;
 
   do{                                // 反復計算のループ
     dx=0.0;
     absx=0.0;
 
     for(i=1;i<=N;i++){
       sum=0;
       for(j=1;j<=N;j++){
 	if(i != j){
 	  sum+=a[i][j]*x[j];
 	}
       }
 
       new=1.0/a[i][i]*(b[i]-sum);   // 反復計算後の近似解
       dx+=fabs(new-x[i]);           // 近似解の変化量を加算
       absx+=fabs(new);              // 近似解の総和計算
       x[i]=new;                     // 新しい近似解を代入
     }
    
   }while(dx/absx > EPS);            // 計算終了条件
  
   for(i=1;i<=N;i++){
     printf("x[%d]=%25.20f\n",i,x[i]);
   }
 
   return 0;
 }

Meta3

Re: 学校の課題について

#7

投稿記事 by Meta3 » 3年前

Windows10     VisualStudio2019   clang

コード:

c:\b>clang-cl c1.c

c:\b>c1.exe
x[1]=   0.99999999999999933387
x[2]=   2.00000000000000044409
x[3]=   3.00000000000000000000

返信

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