C++ オイラー法 わからない

Question

実現したいこと

以下のプログラムを稼働させたいです。

発生している問題・分からないこと

自分はプログラミング初心者で、あまり理解していない部分が多いです。
今回の場合、returnがどのようなものなのかわからず、たぶん意味の分からない使い方をしていると思います。
有識者の方、どのようにしたらよいか教えていただけないでしょうか。

エラーメッセージ

Main.cpp:25:17: warning: unused variable 'x_0' [-Wunused-variable]
         double x_0 = 0.0;//xの初期値
                ^
Main.cpp:59:52: error: called object type 'double' is not a function or function pointer
         result_x = (-v_x/sqrt(v_x*v_x+v_y*v_y))*ρ(v_x*v_x+v_y*v_y)*S*Cd/(2*M) ; 
                                                 ~^
Main.cpp:82:52: error: called object type 'double' is not a function or function pointer
         result_y = (-v_y/sqrt(v_x*v_x+v_y*v_y))*ρ(v_x*v_x+v_y*v_y)*S*Cd/(2*M)-g;
                                                 ~^
1 warning and 2 errors generated.

該当のソースコード

#include <stdio.h> 
 
#include <math.h> 

#define M_PI 3.14159265358979323846
 
double gfunc1(double,double,double,double,double); 

double gfunc2(double,double,double,double,double);
 
double euler2( double, double, double, double, double,double); 
 
  
 
int main(void){ 
 
         double h = 0.01;//オイラー法において設定する刻み幅 
 
         double t_0 =  0.0; //t の初期値 
 
         double v_x_0 =  55*cos((1/2)*M_PI); //Vxの初期値
         
         double v_y_0 =  55*sin((1/2)*M_PI);//Vyの初期値
         
         double x_0 = 0.0;//xの初期値
         
         double y_0 =  0.05; //y の初期値 
 
         double t_n = 2.0; //計算終了時間 
 
  
 
         printf( " result_x = %f [m] \n" , euler2( y_0 , v_x_0 , v_y_0, t_0 , t_n  , h ) ); 
         
         printf( " result_y = %f [m] \n" , euler2( y_0 , v_x_0 ,v_y_0, t_0 , t_n  , h ) ); 
 
         return 0; 
 
} 
 
  
 
double gfunc1( double t, double x , double v_x ,double v_y ){
 
  
 
         double ρ=1.184; 
 
         double M=0.046; 
 
         double S=0.00145;
         
         double Cd=0.35;
         
         double g=9.8066;
 
   double result_x; 
 
         result_x = (-v_x/sqrt(v_x*v_x+v_y*v_y))*ρ(v_x*v_x+v_y*v_y)*S*Cd/(2*M) ; 
 
         return(result_x); 
 
} 

double gfunc2( double t, double x , double v_x ,double v_y ){
 
  
 
         double ρ=1.184; 
 
         double M=0.046; 
 
         double S=0.00145;
         
         double Cd=0.35;
         
         double g=9.8066;
 
   
   double result_y;
         
         result_y = (-v_y/sqrt(v_x*v_x+v_y*v_y))*ρ(v_x*v_x+v_y*v_y)*S*Cd/(2*M)-g;
         
         return(result_y);
 
} 
 
  
 
double euler2( double x_0 , double y_0 ,  double v_x_0, double v_y_0 , double t_0 , double t_n , double h ){ 
 
                         // きざみ幅 h で t_0 での値 y_0,x_0 から t_n まで計算を行い、 y の値を返す関数 
 
  
 
        double t = t_0;//tの初期値を代入する 
        
        double x = y_0;//xの初期値を代入する 
 
        double y = y_0;//yの初期値を代入する 
 
        double v_x = v_x_0;
        
        double v_y = v_y_0;
 
  
        double t_new ;//h 経過後の tの値を保持する 
        
        double x_new ;//h 経過後の xの値を保持する 
 
        double y_new ;//h 経過後の yの値を保持する 
 
        double v_x_new ;
        
        double v_y_new ;
 
  
 
         int i = 0; //ループの回数カウントアップ 
 
  
 
        while( t < t_n ){//tiに対応する変数 t が t_n を超えるまで計算を繰り返す 
 
  
 
         if( i % 10 == 0 ) 
          
                   printf("%f,%f,%f,%f,%f \n" , t , x , y , v_x , v_y  ); 
 
 
  
 
                t_new = t + h; 
 
                v_x_new = v_x + h * gfunc1( t , x , y , v_x, v_y ) ; 
                
                v_y_new = v_y + h * gfunc2( t , x , y , v_x, v_y ) ;
                
                x_new = x + h * v_x; 
 
                y_new = y + h * v_y; 
 
 
  
 
                t = t_new; 
 
                v_x = v_x_new; 
 
                v_y = v_y_new; 
 
                // 次のループで用いることができるよう、値を入れ替える 
 
                   i++; 
 
        } 
 
  
 
        return( x );
        return( y );
 
} 

試したこと・調べたこと

• teratailやGoogle等で検索した
• ソースコードを自分なりに変更した
• 知人に聞いた
• その他

上記の詳細・結果

エラーをコピーして調べたが、よくわからなかった

補足

特になし

Comments

[jimbe]

問題がなければ、なぜこのプログラムを作っているのかを教えていただけますか？
独学/webで c++ を学ばれていて何処かのサイトで見つけた問題たとか、学校/会社から出された課題の穴埋め問題だとか…、

[txty]

return（ｘ）；
return（ｙ）；
とは。一つしか返せないのだから構造体にするとかポインタを引数につかってコピーでなく実体いじればいんじゃないのですかreturnは作業をやめたいとき、関数から外部に取り出したいとき使います。classのメンバにしてもいいのか。

Answer 1

最初のwarningは、x_0が、その後で使用されていないという警告です。
とりあえず、そのままでもよいです。
最終的に使用するか、削除するか自分で決めてください。

残りの2個のエラーは、変数ρの後ろに演算子が無いので、関数として解釈されているというエラーです。
何らかの演算子を入れるか、他の関数に置き換えるかは、数学の問題なので自分で考えてください。

Answer 2

これ、あれだろ。
「空気抵抗がある場合の」斜方投射?を意図してんじゃないか、と。
ρが流体密度、Cdが効力係数、Sが物体の代表面積として、

と言う微分方程式を数値積分(オイラー法)を使って解け、って題意なんじゃないか、って思う。
(チッ、TeratailにはMarkdownはあっても、TeX記法は受け付けないんか・・・)
ホントはその辺、質問文に書いて欲しかったんだが・・・・・・。

んでC++だろ?
僕はC++は専門じゃないんだけど、恐らくC++だとこんな書き方になるんじゃないか、って思う。

#include <iostream>
#include <vector>
#include <tuple>
#include <cmath>
#include <iomanip>
#include <numeric>   // accumulate
#include <algorithm> // reverse
#include <utility>   // pair

using Vec2 = std::pair<double, double>;
using Record = std::tuple<double, Vec2, Vec2>; // (t, r, v)

constexpr double rho = 1.184;
constexpr double Cd  = 0.35;
constexpr double S   = 0.00145;
constexpr double M   = 0.046;
constexpr double g   = 9.8066;
constexpr double h   = 0.01;
constexpr double t_0 = 0.0;
constexpr double t_n = 2.0;

std::vector<double> arange(double start, double end, double step) {
    std::vector<double> result;
    for (double v = start; v < end + 1e-10; v += step)
        result.push_back(v);
    return result;
}

std::vector<Record> foo(const std::vector<Record>& acc, double s) {
    const auto& [t, r, v] = acc.front();  // 先頭が最新状態
    double coeff = - rho * Cd * S * std::sqrt(v.first * v.first + v.second * v.second) / (2 * M);
    Vec2 r_new = {
        r.first + h * v.first,
        r.second + h * v.second
    };
    Vec2 v_new = {
        v.first + h * coeff * v.first,
        v.second + h * (coeff * v.second - g)
    };
    std::vector<Record> result = { std::make_tuple(s, r_new, v_new) };
    result.insert(result.end(), acc.begin(), acc.end());  // prepend
    return result;
}

int main() {
    Vec2 r_0 = {0.0, 0.05};
    Vec2 v_0 = {
        55.0 * std::cos(0.5 * M_PI),
        55.0 * std::sin(0.5 * M_PI)
    };

    std::vector<Record> initial = { std::make_tuple(t_0, r_0, v_0) };
    std::vector<double> steps = arange(t_0 + h, t_n, h);

    std::vector<Record> data = std::accumulate(
        steps.begin(), steps.end(), initial,
        [](const std::vector<Record>& acc, double s) {
            return foo(acc, s);
        }
    );

    std::cout << std::fixed << std::setprecision(10);
    std::reverse(data.begin(), data.end());
    for (const auto& [t, r, v] : data) {
        if (std::fmod(t / 0.01, 10) < 1e-9) {  // t が 0.1, 0.2, 0.3 ... のときだけ出力
            std::cout << t << ", " << r.first << ", " << r.second
                      << ", " << v.first << ", " << v.second << '\n';
        }
    }

    const auto& [tf, rf, vf] = data.back();
    std::cout << " result_x = " << rf.first << " [m]\n";
    std::cout << " result_y = " << rf.second << " [m]\n";

    return 0;
}

「CとC++は似たようなモノ」って誤解があるんだけど、C++にはとても強力なSTLってライブラリが標準添付されていて、それを使ったプログラミングが望ましいスタイルとなるんとちゃうんか。よう知らんけど。
少なくとも知ってる限りで言うと、「C++のエキスパート」はSTLのファンなんだ。だからまずはSTLを使おうとするんじゃないか。異論はあるかもしれんが、ホントその辺はよう知らんのだ。かなり「尖った」C++ファンしか見たことないんで(笑)。
(テンプレートを使ったプログラミングが・・・!!!とか力説されたし、ネットでの某動画の有名なC++な人もそう言ってた気がする)

んで、オイラー法に入る前にまずはプログラミングの基礎的な考え方を。
これは斜方投射なんで二次元空間に於ける記述となる。
こう言った場合、例えばx座標とy座標を「バラバラに扱う」とか言うのってあんま良いテじゃないんだ。
座標情報は座標情報として纏める。
同様に速度も二次元ベクトルなんで、やっぱ要素を「纏める」ってのがセオリーとなるのとちゃうんかな。
ここでは

using Vec2 = std::pair<double, double>;

として二次元ベクトル、Vec2型、ってのを定義している。
厳密に言うと、Vec2型をstd::pair<double, double>へのエイリアスとして用いる、と言う事だ。
これによりr = (x, y)とかv = (vx, vy)みたいな、より物理学らしい表現が可能となる。

もう一つはRecord型の定義だ。

using Record = std::tuple<double, Vec2, Vec2>; // (t, r, v)

これはstd::tuple型へのエイリアスで、含む型はdouble、残りの2つはVec2型と言う前提となっている。
これは貴方のオリジナルのコードだと(時間情報, 位置情報, 速度情報)、と纏めて扱いたかったのかな、ってぇんでこうなってる。コメントに記述された(t, r, v)はそれを意図してる。
そして纏めて扱うってのはモダンなプログラミングに於いては重要な概念なんだ。
もちろん、プログラムをする人の中では、とかく「バラバラにする」とか「分解する」のが好きな人が多い。
多いんだけど、正直言って、そのテの発想は古いプログラミングの考え方なんだよ。
繰り返すけど、多々の情報を纏めて扱うってのはモダンなプログラミングでは重要な概念で、STLを備えたC++はその「モダンなプログラミング」に片足突っ込んでいて、その辺実は古典的なCとは全く違う思想を持っている。
この後それを知ることとなる。

次はユーティリティarangeだ。

std::vector<double> arange(double start, double end, double step) {
    std::vector<double> result;
    for (double v = start; v < end + 1e-10; v += step)
        result.push_back(v);
    return result;
}

これは単純には、この問題上、「時間列」を作るためのユーティリティだ。
貴方の書いたオリジナルのコードでは、恐らく、0〜2秒間を0.01秒刻みで分割したい、んだろう。
言い換えると、0から0.01毎に増えて2に到達する長さ200の数列が欲しい、って事だ。
このユーティリティはそれを満たすstd::vector型の数列を返す関数となる。
なお、ここで言うstd::vectorは上で定義されたVec2型と同じではなく、プログラミング用語では可変長配列を意味してる。

さて、オイラー法に付いて考えていこう。
とは言っても、基本的には、STL上では「考え方」自体は簡単だ。
まずは位置情報rから。
オイラー法では、この問題だと、位置情報rは次の漸化式によって記述される。

ここで、hの単位は実質的には時間(秒)なんで、次元の問題はクリアしてる。
そして速度情報vは

と言う漸化式になっている。
数学的に言うとこの漸化式を解け、ってな問題になるんだろうけど、現象的にはこの、「中途半端な計算過程」に見えるコレ

を、プログラミング用語で畳み込み(folding)と呼ぶ(LaTeXを打ち間違えて無ければいいが・苦笑)。
言い換えると、この「畳み込み」を行う機能をプログラミング言語が提供していれば、この問題は「解ける」って事になるわけだ。C++のSTLではこの「畳み込み」を行うユーティリティ、std::reduceとstd::accumulateが提供されているが、ここではstd::accumulateを使う。

上の「畳み込み」の理論上のモデルだと、なんだかんだ言って「最終値」を把握する事の方が簡単だ。
一方、貴方の書いたオリジナルコードだと、「計算過程」の値も全部見たいらしい。
と言う事は、「計算の過程を全て詰め込んだ情報が必要だ」と言う事になる。
このテの「計算の過程を全て詰め込んだ」データ、あるいはその変数をアキュムレータ(Accumulator: 累積値)と呼ぶ。

今、計算過程の一つが(t, r, v)、つまり「時間、距離、速度」の組だとしてみよう。そしてt_n時のアキュムレータの状態が次のようなモノだと仮定する。

いま、n + 1時点の組を計算する場合、利用すべきはアキュムレータの先頭のデータになる、ってのは自明だろう。そして計算が終わった後、その値はアキュムレータの先頭に追加すべきだ。
これを行うのが関数fooの役割だ。

std::vector<Record> foo(const std::vector<Record>& acc, double s) {
    const auto& [t, r, v] = acc.front();  // 先頭が最新状態
    double coeff = - rho * Cd * S * std::sqrt(v.first * v.first + v.second * v.second);
    Vec2 r_new = {
        r.first + h * v.first,
        r.second + h * v.second
    };
    Vec2 v_new = {
        v.first + h * coeff * v.first,
        v.second + h * (coeff * v.second - g)
    };
    std::vector<Record> result = { std::make_tuple(s, r_new, v_new) };
    result.insert(result.end(), acc.begin(), acc.end());  // prepend
    return result;
}

引数はRecord型vector(可変長引数)のaccとdouble型のsの2つ。
accの先頭からt、r、vの3つの情報を取り出し、そして加速度を出す為の係数(coefficient)を計算する。
そしてaccの先頭情報から位置情報rと速度情報vを計算するわけだが、実はこれは数学的な漸化式をそのまま書けばいいんで何も工夫が必要ない事が分かるだろう。

    Vec2 r_new = {
        r.first + h * v.first,
        r.second + h * v.second
    };
    Vec2 v_new = {
        v.first + h * coeff * v.first,
        v.second + h * (coeff * v.second - g)
    };

あとは新規に(t, r, v)の組を作ってaccの先頭に挿し込む(prepend)してデータ列を返すだけ、だ。
その辺の記述がC++だと若干煩わしいが、「考え方」と言うか、ロジック自体は然程難しくはないだろう。むしろ、この関数だと数式部分が数式記述通りに書けるってメリットは結構デカいんじゃなかろうか。分解しないでいいんだ。

あとはmain関数内でまずは各種初期値を設定する。初期位置、初期速度を設定した後、アキュムレータを変数initialとして定義する。initialはRecord型のstd::vectorだ。
そして時間列stepsを定義する。
アキュムレータと時間列を利用してdataをstd::accumulateで作成する。

    std::vector<Record> data = std::accumulate(
        steps.begin(), steps.end(), initial,
        [](const std::vector<Record>& acc, double s) {
            return foo(acc, s);
        }
    );

時間列がendに到達するまで、std::accumulateは関数fooを利用しつつinitialを加工していく。最終的には時間列の長さ(この問題の場合は200)に合わせて古い順から並んだ200個の計算結果((t, r, v))が詰め込まれたシーケンスを入手する。
これでオイラー法は終了、だ。

あとは出力するだけ、となる。
おしまい。

注1: std::setprecitionではそこに書かれてる通り、「浮動小数点数を出力する精度を設定」してて、std::fixedでは、「浮動小数点数を固定小数点表記で入出力することを指示」しているだけだ。

注2: なお、オリジナルのコードでも速度vのx成分の初期値にcos(1/2 * Π)が使われてるんで、全履歴は水平方向だと位置情報は0のまんま、速度のx成分も0のまんま、となると思う。
結果、事実上このシミュレーションは鉛直方向への「投げ上げ」って事になるんじゃなかろうか。