Top / AA_Experiment / Rocket / SampleProgram_c

AA_Experiment/Rocket/SampleProgram_c

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

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

ソースファイルをダウンロード → filenewton.c

/***********************************************/
/* 航空宇宙工学実験「ロケット推進力の測定」    */
/* 与えられたノズル開口比 Ae/At から、         */
/* 圧力比 pe/pc の数値解をニュートン法で求める */
/***********************************************/

/*
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/pc の初期値 ? ");
    scanf("%lf",&pr);
    printf("\n");
//    pr = 0.1;

    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/pc = %10.8f\n",pr);
}

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

 
Attach file: filenewton.c 95 download [Information]
 
Link: AA_Experiment/Rocket(3469d)
Last-modified: 2005-04-19 (Tue) 12:14:06 (4600d)