1

const int 型の変数がありますが、依存しているパラメーターは double 型です。これを「double」から「const int」にキャストしようとすると、正しく機能しません。たとえば、N が 991 のはずなのに、990 と入力されます。いくつかの方法を試しましたが、1 つしか機能しませんでしたが、この方法が常に機能するかどうかはわかりません。ここに私が試したいくつかの方法があります：

最初の方法:

const int N = (Ls-1)/dx + 1;

2 番目の方法:

const int N = static_cast<const int>((Ls-1)/dx) + 1;

3 番目の方法:

double Z = (Ls-1)/dx + 1;
const int N = Z;

4 番目の方法 (作業方法のみ):

double Z = (Ls-1)/dx;
const int N = Z + 1;

dx は、(Ls-1)/dx の余りが常にゼロになるような値であることに注意してください (つまり、常に整数値です)。型キャストをよりよく理解できるように、他の方法が機能しない理由をとにかく説明できますか?

編集: 要求に応じて、コード全体をアップロードして、すべてがどのように機能しているかを示します。

#include <iostream>
#include <math.h>
#include <stdio.h>
#include <fstream>
#include <cmath>
#include <algorithm>

#define pi 3.14159265

using namespace std;

//Define Fluid Properties
double rho_L = 998; //Liquid Density
double rho_LG = 828.9; //Liquid-Gas Density Ratio
double mu_L = 0.000798; //Liquid Viscosity
double mu_LG = 40.24; //Liquid-Gas Viscosity Ratio
double sigma = 0.0712; //Surface Tension
double nu_G = (mu_L/mu_LG)/(rho_L/rho_LG);

//Define Injector Properties
double Uinj = 56.7; //Injection Velocity
double Dinj = 0.0998; //Injector Diameter
double theta = 15.0*pi/180.0; //Spray Cone Angle
double L = 500.0*Dinj; //Atomization Length
double Ls = L/Dinj; //Normalized Atomization Length

//Define Solver Parameters
double K = 5294; //Viscous Dissipation Coefficient
double Eps = pow(10,-5); //Residual Error
double dx = 0.0001; //Step Size
double Ui = 10; //Initial Guess
//const int Z = static_cast<const int>((Ls-1)/dx + 1) + 1;
const int N = (Ls-1)/dx + 1;//Z;

double deriv (double U, double X, double delta, double m)
{
    double dudx;
    dudx = -(1.0/delta)*(1.0/U)*(U - sqrt(1.0 - U)/sqrt(m*X*X))*(U - sqrt(1.0 - U)/sqrt(m*X*X));
    return (dudx);
}

int main()
{
    //Declare Variables
    int max_step;
    double ERR;
    int step;
    double DEN;
    double SMD;
    double m;
    double Ug;
    double Re;
    double Cd;
    double delta;
    double K1;
    double K2;
    double K3;
    double K4;

    //Allocate Memory From Heap
    double *U = new double [N];
    double *X = new double [N];

    //Initialize Vectors and Variables
    DEN = 0.5*rho_L - (4.0/3.0)*K*(mu_L)/(Uinj*Dinj*Dinj)*L;

    m = 4.0/rho_LG*tan(theta)*tan(theta);

    for (int i = 0; i < N; i++)
    {
        X[i] = 1.0 + dx*i;
    }
    U[0] = 1.0;

    max_step = 1;
    ERR = 1;
    step = 0;
    while(abs(ERR) > Eps && step < max_step)
    {

        //Calculate Ug
        Ug = sqrt(1.0 - (Ui/Uinj))/sqrt(m*Ls*Ls)*Uinj;

        //Calculate SMD
        SMD = 6.0*sigma/(DEN*(Uinj*Uinj - Ui*Ui));

        //Calculate Re # and Drag Coefficient
        Re = abs(Ui-Ug)*SMD/nu_G;

        if(Re <= 0.01)
        {
            Cd = (0.1875) + (24.0/Re);
        }
        else if(Re > 0.01 && Re <= 260.0)
        {
            Cd = (24.0/Re)*(1.0 + 0.1315*pow(Re,(0.32 - 0.05*log10(Re))));
        }
        else
        {
            Cd = (24.0/Re)*(1.0 + 0.1935*pow(Re,0.6305));
        }

        //Determine New U
        delta = (4.0/3.0)*(1.0/Cd)*(rho_LG)*(SMD/Dinj);

        //RK4
        for (int i = 0; i < N-1; i++)
        {
            K1 = deriv(U[i],X[i],delta,m);
            K2 = deriv(U[i]+0.5*dx*K1,X[i]+0.5*dx,delta,m);
            K3 = deriv(U[i]+0.5*dx*K2,X[i]+0.5*dx,delta,m);
            K4 = deriv(U[i]+dx*K3,X[i+1],delta,m);
            U[i+1] = U[i] + dx/6.0*(K1 + 2.0*K2 + 2.0*K3 + K4);
            //if(i >= 0 && i <= 3)
                //cout << i << " " << K1 << " " << K2 << " " << K3 << " " << K4 << " " << U[i] << endl;
        }

        ERR = abs(U[N-1]*Uinj - Ui)/Ui;

        Ui = U[N-1]*Uinj;

        step = step + 1;
    }

    SMD = 6.0*sigma/(DEN*(Uinj*Uinj - Ui*Ui));

    cout << "U = " << Ui << endl;
    cout << "SMD = " << SMD << endl;
    cout << "DEN = " << DEN << endl;
    cout << "Ug = " << Ug << endl;
    cout << "m = " << m << endl;
    cout << "delta = " << delta << endl;
    cout << "Re = " << Re << endl;
    cout << "Cd = " << Cd << endl;
    cout << "U* = " << U[N-1] << endl;
    cout << "Error = " << ERR << endl;
    cout << "step = " << step << endl;

    //Output Data Into Text File
    ofstream outputdata("result-500-15.txt");
    for (int i = 0; i < N; i++)
    {
        outputdata << X[i] << " " << U[i] << '\n';
    }
    outputdata.close();

    delete [] U;
    delete [] X;

    return 0;
}
4

1 に答える 1

3

あなたの推測は正しいです: 0.1 にはバイナリの有限式はありません。これはかなり複雑な問題であり、コメントで述べたように 0.01 を追加しても一般に解決されない多くのコーナーケースがあります。（期待する値などに大きく依存します）

あなたの質問は、商が常に整数であることを示唆しています。doubleその場合、正しい結果を維持するための正しいアプローチは、最初から s を使用しないことです(for Ls、dx、Z)。小数型 (C++ に組み込まれているものは何もなく、独自のライブラリまたはライブラリを使用します)、任意精度の 10 進数型 (ここでも、次のようなライブラリを使用しgmpます。すべての数値に終端の 10 進数式があることがわかっている場合は賢明です)、または、最も簡単: と の両方が小数点以下の桁数であることが保証されている場合は、両方Lsを乗算し、整数型を使用します。dxn10^n

さて、あなたのコードは私が期待したものとは大きく異なります。この場合、私の意見では、正しいことは、ステップ数を修正して、その逆ではなく、そこからN計算することです。dx

const int N = 10000;
double dx = (Ls-1.0)/(double)(N-1);

dx の値から始めて、dx の計算値になるように N を選択する場合は、プログラムの開始時にユーザーに尋ねます。

#include <cmath>

double dxestim;
cout << "dx should be close to: ";
cin >> dxestim;
cout << "Candidate values for N: " << endl;
int N1 = (int) floor((Ls-1)/dx + 1.0);
int N2 = (int) ceil((Ls-1)/dx + 1.0);
cout << N1 << " gives dx = " << (Ls-1.0)/(double)(N1-1) << endl;
cout << N2 << " gives dx = " << (Ls-1.0)/(double)(N2-1) << endl;
cout << "Please choose N: ";
cin >> N;
...
于 2013-02-09T01:09:17.203 に答える