Решение систем диф.уров. методом Рунге-Кутты 4-го порядка

Noob
12 июня 2019
Имеется система из 17-ти дифуров, необходимо решить её методом Рунге-Кутты 4-го порядка.
Я написал код:
    #include <iostream>
#include <cmath>
#include <fstream>
#include <iomanip>


long double f1(long double *k,long double x1, long double x2,long double x3,long double x5,long double x12,long double x13)
{
   return -k[0]*x1 + k[1]*x2*x3 - k[4]*x1*x3 - k[5]*x1*x3 -  k[6]*x1*x5 - k[7]*x1*x5 - k[12]*x1*x2 -  k[13]*x1*x12 - k[16]*x1*x13;
}

long double f2(long double *k,long double x1, long double x2,long double x3,long double x4,long double x5,long double x10,long double x11)
{
    return k[0]*x1 - k[1]*x2*x3 - k[2]*x2 + k[3]*x4*x5 -  k[12]*x1*x2 - k[21]*x10*x2 + k[24]*x11*x5;
}
long double f3(long double *k,long double x1, long double x2,long double x3,long double x4,long double x6,long double x7,long double x9,long double x12,
               long double x13,long double x17)
{
    return k[0]*x1 - k[1]*x2*x3 - k[4]*x1*x3 - k[5]*x1*x3 +  k[8]*x7 - k[17]*x4*x3 + k[18]*x6*x12 - k[19]*x9*x3 +  k[25]*x13 + k[26]*x17 - k[27]*x3*x3 - k[28]*x13*x3 -  k[29]*x12*x3;
}


long double f4(long double *k,long double x1,long double x2, long double x3,long double x4,long double x5,long double x6,long double x7,long double x12,long double x14)
{
    return k[2]*x2 - k[3]*x4*x5 + k[8]*x7 + k[13]*x1*x12 -  k[17]*x4*x3 + k[18]*x6*x12 - k[20]*x4*x12 + k[23]*x14;
}
long double f5(long double *k,long double x1,long double x2,long double x3,long double x4,long double x5,long double x7,long double x8,long double x9,long double x10,
               long double x11,long double x13,long double x14)
{
    return k[2]*x2 - k[3]*x4*x5 - k[6]*x1*x5 - k[7]*x1*x5 +  k[9]*x8 + k[10]*x7 - k[11]*x10*x5 - k[14]*x10*x5 +  k[15]*x9*x13 + k[19]*x9*x3 + k[22]*x14 - k[24]*x11*x5;
}
long double f6(long double *k,long double x1,long double x3,long double x4,long double x6,long double x9,long double x12)
{
    return k[4]*x1*x3 + k[5]*x1*x3 + k[17]*x4*x3 - k[18]*x6*x12 +  k[19]*x9*x3;
}
long double f7(long double *k,long double x1,long double x2,long double x3,long double x5,long double x7,long double x10,long double x13)
{
    return k[4]*x1*x3 + k[6]*x1*x5 - k[8]*x7 - k[10]*x7 + k[11]*x10*x5 + k[12]*x1*x2 + k[16]*x1*x13;
}
long double f8(long double *k,long double x1,long double x3,long double x5,long double x8,long double x12)
{
    return k[5]*x1*x3 + k[7]*x1*x5 - k[9]*x8 + k[13]*x1*x12;
}
long double f9(long double *k,long double x1,long double x3,long double x5,long double x9,long double x10,long double x11,long double x13)
{
    return k[6]*x1*x5 + k[7]*x1*x5 + k[14]*x10*x5 -  k[15]*x9*x13 - k[19]*x9*x3 + k[24]*x11*x5;
}
long double f10(long double *k,long double x1,long double x2,long double x3,long double x5,long double x7,long double x8,long double x9,long double x10,long double x12,long double x13)
{
    return k[9]*x8 + k[10]*x7 - k[11]*x10*x5 - k[14]*x10*x5 +  k[15]*x9*x13 + k[16]*x1*x13 - k[21]*x10*x2 +  k[29]*x12*x3;
}
long double f11(long double *k,long double x1,long double x2,long double x3,long double x5,long double x10,long double x11)
{
    return k[12]*x1*x2 + k[21]*x10*x2 - k[24]*x11*x5 + k[27]*x3*x3;
}
long double f12(long double *k,long double x1,long double x3,long double x4,long double x6,long double x12,long double x14)
{
    return -k[13]*x1*x12 + k[17]*x4*x3 - k[18]*x6*x12 -  k[20]*x4*x12 + k[23]*x14 - k[29]*x12*x3;
}
long double f13(long double *k,long double x1,long double x2,long double x3,long double x5,long double x9,long double x10,long double x13,long double x17)
{
    return k[14]*x10*x5 - k[15]*x9*x13 - k[16]*x1*x13 +  k[21]*x10*x2 - k[25]*x13 + k[26]*x17 - k[28]*x13*x3;
}
long double f14(long double *k,long double x4,long double x12,long double x14)
{
    return k[20]*x4*x12 - k[22]*x14 - k[23]*x14;
}
long double f15(long double *k,long double x14)
{
    return k[22]*x14;
}
long double f16(long double *k,long double x13)
{
    return k[25]*x13;
}
long double f17(long double *k,long double x3,long double x13,long double x17)
{
    return -k[26]*x17 + k[28]*x13*x3;
}

using namespace std;


int main() {
    double h = 0.00001;//шаг
    double t_end=0.0001;
    int n=500;
    float s=0;

    int T=900;
    float R=8.314;
    double *A = new double[30]{ 2.78e10,2.83e5,4.31e4,4.09e4,2.98,5.48,2.55e4,1.13e4,1.20e5,1.60e5,
                                1.09e5,1.30e5,9.70,4.79e2,2.61e3,8.37e2,3.44e2,94.5,128,152,9210,10.2,
                                6.40e4,2.10e5,2.48e4,3e5,1e8,2.64e5,1.64e5,9.56e6};
    double *E = new  double[30]{ 376.,0.,155.,4.15,29.93,22.95,28.27,18.71,126.,150.,149.,13.64,38.25,
                                36.92,10.39,79.49,83.06,39.74,22.86,36.42,
                                19.,27.77,144.2,149.2,35.34,151.,305.,0.,-0.55,0.57};
    long double *k = new  long double[30];
    for(int i=0;i<30;i++){
        k[i]=A[i]*exp(-E[i]/(R*T));
    }
    delete[]A;
    delete[]E;


    long double *x1 = new  long double[n];
    long double *x2 = new  long double[n];
    long  double *x3 = new long  double[n];
    long  double *x4 = new  long double[n];
    long  double *x5 = new  long double[n];
    long  double *x6 = new  long double[n];
    long  double *x7 = new  long double[n];
    long  double *x8 = new  long double[n];
    long  double *x9 = new  long double[n];
    long  double *x10 = new  long double[n];
    long  double *x11 = new  long double[n];
    long  double *x12 = new  long double[n];
    long  double *x13 = new  long double[n];
    long  double *x14 = new  long double[n];
    long  double *x15 = new  long double[n];
    long  double *x16 = new  long double[n];
    long  double *x17 = new  long double[n];


    x1[0]=1,x2[0]=0,x3[0]=0,x4[0]=0,x5[0]=0,x6[0]=0,x7[0]=0,x8[0]=0,x9[0]=0,x10[0]=0,x11[0]=0,x12[0]=0,x13[0]=0,x14[0]=0,x15[0]=0,x16[0]=0,x17[0]=0;
    long double
             /*1*/    a1=0, a2=0, a3=0, a4=0,
             /*2*/    b1=0, b2=0, b3=0, b4=0,
             /*3*/    c1=0, c2=0, c3=0, c4=0,
             /*4*/    d1=0, d2=0, d3=0, d4=0,
             /*5*/    e1=0, e2=0, e3=0, e4=0,
             /*6*/    F1=0, F2=0, F3=0, F4=0,
             /*7*/    g1=0, g2=0, g3=0, g4=0,
             /*8*/    h1=0, h2=0, h3=0, h4=0,
             /*9*/    i1=0, i2=0, i3=0, i4=0,
             /*10*/   j1=0, j2=0, j3=0, j4=0,
             /*11*/   k1=0, k2=0, k3=0, k4=0,
             /*12*/   l1=0, l2=0, l3=0, l4=0,
             /*13*/   m1=0, m2=0, m3=0, m4=0,
             /*14*/   n1=0, n2=0, n3=0, n4=0,
             /*15*/   o1=0, o2=0, o3=0, o4=0,
             /*16*/   p1=0, p2=0, p3=0, p4=0,
             /*17*/   q1=0, q2=0, q3=0, q4=0;
    int i=1;

    for(double s=0;s<t_end;s+=h){
        a1 = h*f1(k,x1[i-1],x2[i-1],x3[i-1],x5[i-1],x12[i-1],x13[i-1]);
        b1 = h*f2(k,x1[i-1],x2[i-1],x3[i-1],x4[i-1],x5[i-1],x10[i-1],x11[i-1]);
        c1 = h*f3(k,x1[i-1],x2[i-1],x3[i-1],x4[i-1],x6[i-1],x7[i-1],x9[i-1],x12[i-1],x13[i-1],x17[i-1]);
        d1 = h*f4(k,x1[i-1],x2[i-1],x3[i-1],x4[i-1],x5[i-1],x6[i-1],x7[i-1],x12[i-1],x14[i-1]);
        e1 = h*f5(k,x1[i-1],x2[i-1],x3[i-1],x4[i-1],x5[i-1],x7[i-1],x8[i-1],x9[i-1],x10[i-1],x11[i-1],x13[i-1],x14[i-1]);
        F1 = h*f6(k,x1[i-1],x3[i-1],x4[i-1],x6[i-1],x9[i-1],x12[i-1]);
        g1 = h*f7(k,x1[i-1],x2[i-1],x3[i-1],x5[i-1],x7[i-1],x10[i-1],x13[i-1]);
        h1 = h*f8(k,x1[i-1],x3[i-1],x5[i-1],x8[i-1],x12[i-1]);
        i1 = h*f9(k,x1[i-1],x3[i-1],x5[i-1],x9[i-1],x10[i-1],x11[i-1],x13[i-1]);
        j1 = h*f10(k,x1[i-1],x2[i-1],x3[i-1],x5[i-1],x7[i-1],x8[i-1],x9[i-1],x10[i-1],x12[i-1],x13[i-1]);
        k1 = h*f11(k,x1[i-1],x2[i-1],x3[i-1],x5[i-1],x10[i-1],x11[i-1]);
        l1 = h*f12(k,x1[i-1],x3[i-1],x4[i-1],x6[i-1],x12[i-1],x14[i-1]);
        m1 = h*f13(k,x1[i-1],x2[i-1],x3[i-1],x5[i-1],x9[i-1],x10[i-1],x13[i-1],x17[i-1]);
        n1 = h*f14(k,x4[i-1],x12[i-1],x14[i-1]);
        o1 = h*f15(k,x14[i-1]);
        p1 = h*f16(k,x13[i-1]);
        q1 = h*f17(k,x3[i-1],x13[i-1],x17[i-1]);




        a2 = h*f1(k,x1[i-1]+a1/2,x2[i-1]+b1/2,x3[i-1]+c1/2,x5[i-1]+e1/2,x12[i-1]+l1/2,x13[i-1]+m1/2);
        b2 = h*f2(k,x1[i-1]+a1/2,x2[i-1]+b1/2,x3[i-1]+c1/2,x4[i-1]+d1/2,x5[i-1]+e1/2,x10[i-1]+j1/2,x11[i-1]+k1/2);
        c2 = h*f3(k,x1[i-1]+a1/2,x2[i-1]+b1/2,x3[i-1]+c1/2,x4[i-1]+d1/2,x6[i-1]+F1/2,x7[i-1]+g1/2,x9[i-1]+i1/2,x12[i-1]+l1/2,x13[i-1]+m1/2,x17[i-1]+q1/2);
        d2 = h*f4(k,x1[i-1]+a1/2,x2[i-1]+b1/2,x3[i-1]+c1/2,x4[i-1]+d1/2,x5[i-1]+e1/2,x6[i-1]+F1/2,x7[i-1]+g1/2,x12[i-1]+l1/2,x14[i-1]+n1/2);
        e2 = h*f5(k,x1[i-1]+a1/2,x2[i-1]+b1/2,x3[i-1]+c1/2,x4[i-1]+d1/2,x5[i-1]+e1/2,x7[i-1]+g1/2,x8[i-1]+h1/2,x9[i-1]+i1/2,x10[i-1]+j1/2,x11[i-1]+k1/2,x13[i-1]+m1/2,x14[i-1]+n1/2);
        F2 = h*f6(k,x1[i-1]+a1/2,x3[i-1]+c1/2,x4[i-1]+d1/2,x6[i-1]+F1/2,x9[i-1]+i1/2,x12[i-1]+l1/2);
        g2 = h*f7(k,x1[i-1]+a1/2,x2[i-1]+b1/2,x3[i-1]+c1/2,x5[i-1]+e1/2,x7[i-1]+g1/2,x10[i-1]+j1/2,x13[i-1]+m1/2);
        h2 = h*f8(k,x1[i-1]+a1/2,x3[i-1]+c1/2,x5[i-1]+e1/2,x8[i-1]+h1/2,x12[i-1]+l1/2);
        i2 = h*f9(k,x1[i-1]+a1/2,x3[i-1]+c1/2,x5[i-1]+e1/2,x9[i-1]+i1/2,x10[i-1]+j1/2,x11[i-1]+k1/2,x13[i-1]+m1/2);
        j2 = h*f10(k,x1[i-1]+a1/2,x2[i-1]+b1/2,x3[i-1]+c1/2,x5[i-1]+e1/2,x7[i-1]+g1/2,x8[i-1]+h1/2,x9[i-1]+i1/2,x10[i-1]+j1/2,x12[i-1]+l1/2,x13[i-1]+m1/2);
        k2 = h*f11(k,x1[i-1]+a1/2,x2[i-1]+b1/2,x3[i-1]+c1/2,x5[i-1]+e1/2,x10[i-1]+j1/2,x11[i-1]+k1/2);
        l2 = h*f12(k,x1[i-1]+a1/2,x3[i-1]+c1/2,x4[i-1]+d1/2,x6[i-1]+F1/2,x12[i-1]+l1/2,x14[i-1]+n1/2);
        m2 = h*f13(k,x1[i-1]+a1/2,x2[i-1]+b1/2,x3[i-1]+c1/2,x5[i-1]+e1/2,x9[i-1]+i1/2,x10[i-1]+j1/2,x13[i-1]+m1/2,x17[i-1]+q1/2);
        n2 = h*f14(k,x4[i-1]+d1/2,x12[i-1]+l1/2,x14[i-1]+n1/2);
        o2 = h*f15(k,x14[i-1]+n1/2);
        p2 = h*f16(k,x13[i-1]+m1/2);
        q2 = h*f17(k,x3[i-1]+c1/2,x13[i-1]+m1/2,x17[i-1]+q1/2);


        a3 = h*f1(k,x1[i-1]+a2/2,x2[i-1]+b2/2,x3[i-1]+c2/2,x5[i-1]+e2/2,x12[i-1]+l2/2,x13[i-1]+m2/2);
        b3 = h*f2(k,x1[i-1]+a2/2,x2[i-1]+b2/2,x3[i-1]+c2/2,x4[i-1]+d2/2,x5[i-1]+e2/2,x10[i-1]+j2/2,x11[i-1]+k2/2);
        c3 = h*f3(k,x1[i-1]+a2/2,x2[i-1]+b2/2,x3[i-1]+c2/2,x4[i-1]+d2/2,x6[i-1]+F2/2,x7[i-1]+g2/2,x9[i-1]+i2/2,x12[i-1]+l2/2,x13[i-1]+m2/2,x17[i-1]+q2/2);
        d3 = h*f4(k,x1[i-1]+a2/2,x2[i-1]+b2/2,x3[i-1]+c2/2,x4[i-1]+d2/2,x5[i-1]+e2/2,x6[i-1]+F2/2,x7[i-1]+g2/2,x12[i-1]+l2/2,x14[i-1]+n2/2);
        e3 = h*f5(k,x1[i-1]+a2/2,x2[i-1]+b2/2,x3[i-1]+c2/2,x4[i-1]+d2/2,x5[i-1]+e2/2,x7[i-1]+g2/2,x8[i-1]+h2/2,x9[i-1]+i2/2,x10[i-1]+j2/2,x11[i-1]+k2/2,x13[i-1]+m2/2,x14[i-1]+n2/2);
        F3 = h*f6(k,x1[i-1]+a2/2,x3[i-1]+c2/2,x4[i-1]+d2/2,x6[i-1]+F2/2,x9[i-1]+i2/2,x12[i-1]+l2/2);
        g3 = h*f7(k,x1[i-1]+a2/2,x2[i-1]+b2/2,x3[i-1]+c2/2,x5[i-1]+e2/2,x7[i-1]+g2/2,x10[i-1]+j2/2,x13[i-1]+m2/2);
        h3 = h*f8(k,x1[i-1]+a2/2,x3[i-1]+c2/2,x5[i-1]+e2/2,x8[i-1]+h2/2,x12[i-1]+l2/2);
        i3 = h*f9(k,x1[i-1]+a2/2,x3[i-1]+c2/2,x5[i-1]+e2/2,x9[i-1]+i2/2,x10[i-1]+j2/2,x11[i-1]+k2/2,x13[i-1]+m2/2);
        j3 = h*f10(k,x1[i-1]+a2/2,x2[i-1]+b2/2,x3[i-1]+c2/2,x5[i-1]+e2/2,x7[i-1]+g2/2,x8[i-1]+h2/2,x9[i-1]+i2/2,x10[i-1]+j2/2,x12[i-1]+l2/2,x13[i-1]+m2/2);
        k3 = h*f11(k,x1[i-1]+a2/2,x2[i-1]+b2/2,x3[i-1]+c2/2,x5[i-1]+e2/2,x10[i-1]+j2/2,x11[i-1]+k2/2);
        l3 = h*f12(k,x1[i-1]+a2/2,x3[i-1]+c2/2,x4[i-1]+d2/2,x6[i-1]+F2/2,x12[i-1]+l2/2,x14[i-1]+n2/2);
        m3 = h*f13(k,x1[i-1]+a2/2,x2[i-1]+b2/2,x3[i-1]+c2/2,x5[i-1]+e2/2,x9[i-1]+i2/2,x10[i-1]+j2/2,x13[i-1]+m2/2,x17[i-1]+q2/2);
        n3 = h*f14(k,x4[i-1]+d2/2,x12[i-1]+l2/2,x14[i-1]+n2/2);
        o3 = h*f15(k,x14[i-1]+n2/2);
        p3 = h*f16(k,x13[i-1]+m2/2);
        q3 = h*f17(k,x3[i-1]+c2/2,x13[i-1]+m2/2,x17[i-1]+q2/2);


        a4 = h*f1(k,x1[i-1]+a3,x2[i-1]+b3,x3[i-1]+c3,x5[i-1]+e3,x12[i-1]+l3,x13[i-1]+m3);
        b4 = h*f2(k,x1[i-1]+a3,x2[i-1]+b3,x3[i-1]+c3,x4[i-1]+d3,x5[i-1]+e3,x10[i-1]+j3,x11[i-1]+k3);
        c4 = h*f3(k,x1[i-1]+a3,x2[i-1]+b3,x3[i-1]+c3,x4[i-1]+d3,x6[i-1]+F3,x7[i-1]+g3,x9[i-1]+i3,x12[i-1]+l3,x13[i-1]+m3,x17[i-1]+q3);
        d4 = h*f4(k,x1[i-1]+a3,x2[i-1]+b3,x3[i-1]+c3,x4[i-1]+d3,x5[i-1]+e3,x6[i-1]+F3,x7[i-1]+g3,x12[i-1]+l3,x14[i-1]+n3);
        e4 = h*f5(k,x1[i-1]+a3,x2[i-1]+b3,x3[i-1]+c3,x4[i-1]+d3,x5[i-1]+e3,x7[i-1]+g3,x8[i-1]+h3,x9[i-1]+i3,x10[i-1]+j3,x11[i-1]+k3,x13[i-1]+m3,x14[i-1]+n3);
        F4 = h*f6(k,x1[i-1]+a3,x3[i-1]+c3,x4[i-1]+d3,x6[i-1]+F3,x9[i-1]+i3,x12[i-1]+l3);
        g4 = h*f7(k,x1[i-1]+a3,x2[i-1]+b3,x3[i-1]+c3,x5[i-1]+e3,x7[i-1]+g3,x10[i-1]+j3,x13[i-1]+m3);
        h4 = h*f8(k,x1[i-1]+a3,x3[i-1]+c3,x5[i-1]+e3,x8[i-1]+h3,x12[i-1]+l3);
        i4 = h*f9(k,x1[i-1]+a3,x3[i-1]+c3,x5[i-1]+e3,x9[i-1]+i3,x10[i-1]+j3,x11[i-1]+k3,x13[i-1]+m3);
        j4 = h*f10(k,x1[i-1]+a3,x2[i-1]+b3,x3[i-1]+c3,x5[i-1]+e3,x7[i-1]+g3,x8[i-1]+h3,x9[i-1]+i3,x10[i-1]+j3,x12[i-1]+l3,x13[i-1]+m3);
        k4 = h*f11(k,x1[i-1]+a3,x2[i-1]+b3,x3[i-1]+c3,x5[i-1]+e3,x10[i-1]+j3,x11[i-1]+k3);
        l4 = h*f12(k,x1[i-1]+a3,x3[i-1]+c3,x4[i-1]+d3,x6[i-1]+F3,x12[i-1]+l3,x14[i-1]+n3);
        m4 = h*f13(k,x1[i-1]+a3,x2[i-1]+b3,x3[i-1]+c3,x5[i-1]+e3,x9[i-1]+i3,x10[i-1]+j3,x13[i-1]+m3,x17[i-1]+q3);
        n4 = h*f14(k,x4[i-1]+d3,x12[i-1]+l3,x14[i-1]+n3);
        o4 = h*f15(k,x14[i-1]+n3);
        p4 = h*f16(k,x13[i-1]+m3);
        q4 = h*f17(k,x3[i-1]+c3,x13[i-1]+m3,x17[i-1]+q3);

        x1[i] =  x1[i-1]+(a1+2*a2+2*a3+a4)/6;
        x2[i] =  x2[i-1]+(b1+2*b2+2*b3+b4)/6;
        x3[i] =  x3[i-1]+(c1+2*c2+2*c3+c4)/6;
        x4[i] =  x4[i-1]+(d1+2*d2+2*d3+d4)/6;
        x5[i] =  x5[i-1]+(e1+2*e2+2*e3+e4)/6;
        x6[i] =  x6[i-1]+(F1+2*F2+2*F3+F4)/6;
        x7[i] =  x7[i-1]+(g1+2*g2+2*g3+g4)/6;
        x8[i] =  x8[i-1]+(h1+2*h2+2*h3+h4)/6;
        x9[i] =  x9[i-1]+(i1+2*i2+2*i3+i4)/6;
        x10[i] = x10[i-1]+(j1+2*j2+2*j3+j4)/6;
        x11[i] = x11[i-1]+(k1+2*k2+2*k3+k4)/6;
        x12[i] = x12[i-1]+(l1+2*l2+2*l3+l4)/6;
        x13[i] = x13[i-1]+(m1+2*m2+2*m3+m4)/6;
        x14[i] = x14[i-1]+(n1+2*n2+2*n3+n4)/6;
        x15[i] = x15[i-1]+(o1+2*o2+2*o3+o4)/6;
        x16[i] = x16[i-1]+(p1+2*p2+2*p3+p4)/6;
        x17[i] = x17[i-1]+(q1+2*q2+2*q3+q4)/6;
        i++;

    }
    cout<<" t        x1"<<endl;
    double q=0;
  for(int j=0;j<i;j++){
    cout<<q<<" "<<x1[j]<<endl;
     q+=h;
    }
    delete[] x1;
    delete[] x2;
    delete[] x3;
    delete[] x4;
    delete[] x5;
    delete[] x6;
    delete[] x7;
    delete[] x8;
    delete[] x9;
    delete[] x10;
    delete[] x11;
    delete[] x12;
    delete[] x13;
    delete[] x14;
    delete[] x15;
    delete[] x16;
    delete[] x17;



    return 0;
}
Но к сожалению, он не работает. Начальные данные верные. Проверял несколько раз, ошибка скорее в алгоритме.
Должны получиться такие значения x1:
{1., 2.0360510^-6, 9.6865610^-7, 6.3547310^-7, 4.8322510^-7, 3.949710^-7, 3.3594810^-7, 2.9300310^-7, 2.6007110^-7, 2.3390410^-7, 2.1256310^-7}
Решение систем диф.уров. методом Рунге-Кутты 4-го порядка

Ответить
