В ответе inf и nan (решений уравнений методом Эйлера).

[KEIego]

Здравствуйте!

При реализации задачи на С++, в ответе получила inf (как я понимаю, бесконечность) и nan. Возможно, влияет то, что нужно для начала инициализировать все переменные (включая Ucbf, Uc1f...)? Или необходимо выбирать более мелкий шаг (dt)? Решение легко нашлось в системе Wolfram Mathematica, по ним строились графики.
Помогите, пожалуйста, найти проблему!
Думаю, что метод Эйлера реализован правильно (решение проверялось на закомментированной системе //Ucf = Uci + Ili/C*(tf-ti); //Ilf = Ili + (E-Uci)/L*(tf-ti)

Исходные данные:

E2 = 5;
Cb = 2 10^-12;
C1 = 10^-6;
C2 = 10^-6;
i = 10^-12;
R = 10^6;
Rb = 20;
R1 = 1000;
R2 = 1000;
R3 = 1000;
L1 = 10^-5;
E1[t] := 10 Sin[20000 Pi t];

Система:

UCb'[t] == (1/Cb) ((1/Rb) (E1[t] - E2 - UC1[t] - UCb[t]) - i (Exp[UCb[t]/0.026 - 1]) - (1/R) UCb[t]),
UC1'[t] == (1/C1) (IL1[t] + (1/Rb) (E1[t] - E2 - UC1[t] - UCb[t]) - (1/R1) UC1[t]),
UC2'[t] == (1/C2) (IL1[t] - (1/R3) UC2[t]),
IL1'[t] == (-1/L1) (UC1[t] + UC2[t] - E1[t] + R2 IL1[t]),
UCb[0] == 0,
UC1[0] == 0,
UC2[0] == 0,
IL1[0] == 0.

#include <iostream>
#include <fstream>
#include <iomanip>
#include <cmath>
#include <string>
 
#define W_E1 3.1415926535897932384626433832795
#define E2 5.0
#define Cb 2e-12
#define C1 1e-6
#define C2 1e-6
#define I 1e-12
#define R 1e+6
#define Rb 20.0
#define R1 1000.0
#define R2 1000.0
#define R3 1000.0
#define L 1e-5
 
using namespace std;
 
double euler(double, double, double, double, double, double, double, double&, double&, double&, double&); //передача параметров по ссылке
 
int main()
{
    double E1i,E10,ti, tf, dt, tmax, Ucbf, Uc1f, Uc2f, Ilf, Ucbi, Uc1i, Uc2i, Ili;
    const string method = "Simple Euler";
 
/* output: file and formats */
    ofstream file;
    file.open ("ode.txt");
    file.precision(6);
    file.setf(ios::fixed | ios::showpoint);
    cout.precision(6);
    cout.setf(ios::fixed | ios::showpoint);
 
/* initial information */
    ti = 0.0;             // initial value for variable
    dt = 0.00001;             // step size for integration
    tmax = 0.001 ;          // inegrate from ti till tmax
    Ucbi = 0.0;
    Uc1i = 0.0;
    Uc2i = 0.0;
    Ili = 0.0;
    E10 = 10*sin(20000*W_E1*ti);
 
 
    file << setw(24) << method << endl;
    file << setw(16) << "t" << setw(16) << "ucb" << setw(16) << "uc1" << setw(16) << "uc2" << setw(16) << "il" << setw(16) << "E1" << endl;
    file << setw(16) << ti << setw(16) << Ucbi << setw(16) << Uc1i << setw(16) << Uc2i << setw(16) << Ili << setw(16) << E10 << endl;
 
/* step 2: integration of ODE */
    while (ti <= tmax)
    {
        tf = ti + dt;
        E1i = 10*sin(20000*W_E1*(ti+dt));
        euler(E1i,ti,tf,Ucbi,Uc1i,Uc2i,Ili,Ucbf,Uc1f,Uc2f,Ilf);
        file << setw(16) << tf << setw(16) << Ucbf << setw(16) << Uc1f << setw(16) << Uc2f << setw(16) << Ilf << setw(16) << E1i << endl;
        ti = tf;
        Ucbi = Ucbf;
        Uc1i = Uc1f;
        Uc2i = Uc2f;
        Ili = Ilf;
    }
 
    cout << "Hello World!" << endl;
    return 0;
}
 
double euler(double E1i, double ti, double tf, double Ucbi, double Uc1i, double Uc2i, double Ili, double& Ucbf, double& Uc1f, double& Uc2f, double& Ilf)
{
     //Ucf = Uci + Ili/C*(tf-ti);
     //Ilf = Ili + (E-Uci)/L*(tf-ti);
    Ucbf = Ucbi + ((1/Rb)*(E1i - E2 - Uc1i - Ucbi) - (1/R)*Ucbi - I*exp(Ucbi/0.026-1))/Cb*(tf-ti);
    Uc1f = Uc1i + (Ili + (1/Rb)*(E1i - E2 - Uc1i - Ucbi) - (1/R1)*Uc1i)/C1*(tf-ti);
    Uc2f = Uc2i + (Ili - (1/R3)*Uc2i)/C2*(tf-ti);
    Ilf = Ili - (Uc1i + Uc2i - E1i + R2*Ili)/L*(tf-ti);
}

Вывод значений в файле:
t ucb uc1 uc2 il E1
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
0.000010 219463.130731 0.438926 0.000000 5.877853 5.877853
0.000020 -inf -109670.316484 58.778525 -5862.903032 9.510565
0.000030 nan inf -58570.839575 5966661.176993 9.510565
0.000040 nan nan 59608626.638747 -inf 5.877853
0.000050 nan nan -inf nan -0.000000
0.000060 nan nan nan nan -5.877853
0.000070 nan nan nan nan -9.510565
0.000080 nan nan nan nan -9.510565
0.000090 nan nan nan nan -5.877853
0.000100 nan nan nan nan 0.000000
И т.д.