BOOST:ODEINT突然迭代停止

问题描述

我是C ++领域的新手,我在boost库上遇到了一些麻烦。在我的问题中,我想用5个方程式求解ODE系统。这不是一个僵化的问题。作为迭代方法,我同时使用了integreate(rhs,x0,t0,tf,size_step,write_output)integreate_adaptive(stepper,sys,write_output)。这两种方法实际上都整合了方程,但是给了我无意义的结果,几乎将步长从0.001随机更改为5。方程和数据正确。我该怎么做才能解决此问题?这是代码

#include <iostream>
#include <vector>
#include <boost/numeric/odeint.hpp>
#include <fstream>
#include <boost/array.hpp>

using namespace std;
using namespace boost::numeric::odeint;

//DATA
double Lin = 20000;                             // kg/h
double Gdry = 15000;                            // kg/h
double P = 760;                                 // mmHg
double TinH2O = 50;                             // °C
double ToutH2O = 25;                            // °C
double Tinair = 20;                             // °C
double Z = 0.5;                                 // relative humidity
double Cu = 0.26;                               // kcal/kg*K
double CpL = 1;                                 // kcal/kg*K
double DHev = 580;                              // kcal/kg
double hga = 4000;                              // kcal/m/h/K
double hla = 30000;                             // kcal/m/h/K
double A = -49.705;                             // Pev 1st coeff    mmHg vs °C
double B = 2.71;                                // Pev 2nd coeff    mmHg vs °C
double Usair = 0.62*(A + B*Tinair) / P;
double Uair = Z*Usair;
double Kua = hga / Cu;
double L0 = 19292;                              // kg/h


typedef vector< double > state_type;
vector <double> pack_height;
vector <double> Umidity;
vector <double> T_liquid;
vector <double> T_gas;
vector <double> Liquid_flow;
vector <double> Gas_flow;

void rhs(const state_type& x,state_type& dxdt,const double z )
{// U Tl Tg L G

double Ti = (hla*x[1] + hga*x[2] + Kua*DHev*(x[0] - 0.62*A / P)) / (hla + hga + Kua*DHev*0.62*B / P);
double Ui = 0.62*(A + B*Ti) / P;

dxdt[0] = Kua*(Ui - x[0]) / Gdry / 100;
dxdt[1] = hla*(x[1] - Ti) / x[3] / CpL / 100;
dxdt[2] = hga*(Ti - x[2]) / Gdry / Cu / 100;
dxdt[3] = Kua*(Ui - x[0]) / 100;
dxdt[4] = Kua*(Ui - x[0]) / 100;
}

void write_output(const state_type& x,const double z)
{
pack_height.push_back(z);
Umidity.push_back(x[0]);
T_liquid.push_back(x[1]);
T_gas.push_back(x[2]);
Liquid_flow.push_back(x[3]);
Gas_flow.push_back(x[4]);

cout << z << "      " << x[0] << "      " << x[1] << "      " << x[2] << "      " << x[3] << "      " << x[4] << endl;
}

int main()
{
state_type x(5);
x[0] = Uair;
x[1] = ToutH2O;
x[2] = Tinair;
x[3] = L0;
x[4] = Gdry;

double z0 = 0.0;
double zf = 5.5;
double stepsize = 0.001;
integrate( rhs,x,z0,zf,stepsize,write_output );

return 0;
}

这是我从提示符处得到的最终结果:

    0         0.00183349      25           20           19292      15000
    0.001     0.00183356      25           20           19292      15000
    0.0055    0.0018339       25.0002      20.0001      19292      15000
    0.02575   0.00183542      25.001       20.0007      19292      15000
    0.116875  0.00184228      25.0046      20.003       19292.1    15000.1
    0.526938  0.00187312      25.0206      20.0135      19292.6    15000.6
    2.37222   0.00201203      25.0928      20.0608      19294.7    15002.7
    5.5       0.00224788      25.2155      20.142       19298.2    15006.2

只有第一个迭代具有正确的步进大小..很明显,解决方案不是正确的..我该怎么办?先感谢您。 :)

解决方法

如果您read the documentation,则将发现恒定步长例程为integrate_constintegrate_n_steps,或者可能是integrate_adaptive,且步长不受控制。 / p>

integrate的简短调用使用具有自适应步长的标准dopri5步进器,因此步长的变化并不奇怪。您可以使用步进器的密集输出在等距时间内插值。