Top / AA_Experiment / Rocket / SampleProgram_c

Diff of AA_Experiment/Rocket/SampleProgram_c


*圧力比 pe/p0 の数値解をニュートン法で求めるサンプルプログラム (C言語版) [#z4fa9462]
*圧力比 pe/pc の数値解をニュートン法で求めるサンプルプログラム (C言語版) [#z4fa9462]

下記の内容で、拡張子".c"のファイルを作成し、Cコンパイラでコンパイルし、実行する。

ソースファイルをダウンロード → &ref(newton.c);

 /***********************************************/
 /* 航空宇宙工学実験「ロケット推進力の測定」    */
 /* 与えられたノズル開口比 Ae/At から、         */
 /* 圧力比 pe/p0 の数値解をニュートン法で求める */
 /* 圧力比 pe/pc の数値解をニュートン法で求める */
 /***********************************************/
 
 /*
 pe/p0=x とおき、問題の式を整理して、
 pe/pc=x とおき、問題の式を整理して、
 
         2/γ   (γ+1)/γ         2                        (γ+1)/(γ-1)
 f(x) = x    - x         - (At/Ae) * (γ-1)/2) * (2/(γ+1))              = 0
 
 を満たす x を求める。
 */
 
 
 #include <stdio.h>
 #include <math.h>
 
 /*  定数  */
 #define PI  3.14159265 //円周率π
 #define EPS 1.0e-10    //収束条件ε
 #define G   1.402      //比熱比γ
 
 double fx(double x , double ar) {
     return pow(x,(2.0/G))-pow(x,(1.0+1.0/G))
         -((G-1.0)/2.0*pow(2.0/(G+1.0),(G+1.0)/(G-1.0)))/ar/ar;
 }
 
 double dfx(double x) {
     return (2.0/G) * pow(x,(2.0/G-1.0))
         -(1.0+1.0/G) * pow(x,(1.0/G));
 }
 
 void main(void) {
     int i=0;
     double dt, de, at, ae, ar, pr, f, df, dx=1.0;
 
     /* ノズル寸法の入力 */
     printf("スロート直径 Dt ? ");
     scanf("%lf",&dt);
     printf("  出口直径 De ? ");
     scanf("%lf",&de);
 //    dt = 2.5;
 //    de = 3.3;
     at = dt*dt*PI/4.0;
     ae = de*de*PI/4.0;
     ar = (de*de/dt/dt);
     printf("\n");
     printf("スロート面積 At = %g\n", at);
     printf("  出口面積 Ae = %g\n", ae);
     printf("開口比  Ae/At = %g\n\n", ar);
 
     /* ニュートン法を開始する初期値 */
     printf("圧力比 pe/p0 の初期値 ? ");
     printf("圧力比 pe/pc の初期値 ? ");
     scanf("%lf",&pr);
     printf("\n");
 //    pr = 0.1;
 
     printf(" i   pe/p0        f              df\n");
     printf(" i   pe/pc        f              df\n");
     while(fabs(dx) > EPS) {    // 収束するまで繰り返し
         f = fx(pr,ar);
         df = dfx(pr);
         printf("%2d  %10.8f  %10g  %10g\n",i,pr,f,df);
         dx = - f/df;
         pr += dx;
         i++;
     }
     printf("\n");
 
     /* 収束値の表示 */
     printf("圧力比 pe/p0 = %10.8f\n",pr);
     printf("圧力比 pe/pc = %10.8f\n",pr);
 }

理解しやすさを優先して、かなり無駄な書き方をしています。