IIshizaka and Flanagan's two-mass model of the vocal cords.

2mass.cpp

#include <iostream>
#include <cstdlib>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_odeiv.h>
#include <cmath>


const gsl_odeiv_step_type *StepType = gsl_odeiv_step_rk4;
const double dt = 1.E-5;
const double Hmin = 1.E-3;
const double T = 1.0;

const double m1  = 0.125;
const double m2  = 0.025;
const double ks1 = 80000;
const double ks2 = 8000;
const double es1 = 100;
const double es2 = 100;
const double kh1 = 3 * ks1;
const double kh2 = 3 * ks2;
const double eh1 = 500;
const double eh2 = 500;
const double kc  = 25000;
const double om1 = std::sqrt(ks1 / m1);
const double om2 = std::sqrt(ks2 / m2);
const double xi1 = 0.1;
const double xi2 = 0.6;
const double lg  = 1.4;
const double d1  = 0.25;
const double d2  = 0.05;
const double Ag0 = 0.05;
const double Ps  = 7850;
const double xc = -Ag0 / (2 * lg);

//constexpr double A[] = { 0.45, 0.20, 0.26, 0.21, 0.32, 0.30, 0.33, 1.05, 1.12, 0.85, 0.63, 0.39, 0.26, 0.28, 0.23, 0.32, 0.29, 0.28, 0.40, 0.66, 1.20, 1.05, 1.62, 2.09, 2.56, 2.78, 2.86, 3.02, 3.75, 4.60, 5.09, 6.02, 6.55, 6.29, 6.27, 5.94, 5.28, 4.70, 3.87, 4.13, 4.25, 4.27, 4.69, 5.03 }; const char tract[] = "a";
//constexpr double A[] = { 0.33, 0.30, 0.36, 0.34, 0.68, 0.50, 2.43, 3.15, 2.66, 2.49, 3.39, 3.80, 3.78, 4.35, 4.50, 4.43, 4.68, 4.52, 4.15, 4.09, 3.51, 2.95, 2.03, 1.66, 1.38, 1.05, 0.60, 0.35, 0.32, 0.12, 0.10, 0.16, 0.25, 0.24, 0.38, 0.28, 0.36, 0.65, 1.58, 2.05, 2.01, 1.58 }; const char[] tract = "i";
//constexpr double A[] = { 0.40, 0.38, 0.28, 0.43, 0.55, 1.72, 2.91, 2.88, 2.37, 2.10, 3.63, 5.86, 5.63, 5.43, 4.80, 4.56, 4.29, 3.63, 3.37, 3.16, 3.31, 3.22, 2.33, 2.07, 2.07, 1.52, 0.74, 0.23, 0.15, 0.22, 0.22, 0.37, 0.60, 0.76, 0.86, 1.82, 2.35, 2.55, 3.73, 5.47, 4.46, 2.39, 1.10, 0.77, 0.41, 0.86 }; const char tract[] = "u";
//constexpr double A[] = { 0.21, 0.13, 0.16, 0.14, 0.06, 0.78, 1.25, 1.24, 0.99, 0.72, 0.73, 1.06, 1.77, 1.97, 2.46, 2.70, 2.92, 3.03, 2.84, 2.84, 2.83, 2.36, 2.14, 2.00, 1.78, 1.81, 1.79, 1.50, 1.37, 1.36, 1.43, 1.83, 2.08, 2.59, 2.54, 2.11, 2.34, 2.74, 2.19, 1.60 }; const char tract[] = "e";
constexpr double A[] = { 0.18, 0.17, 0.23, 0.28, 0.59, 1.46, 1.60, 1.11, 0.82, 1.01, 2.72, 2.71, 1.96, 1.92, 1.70, 1.66, 1.52, 1.28, 1.44, 1.28, 0.89, 1.25, 1.38, 1.09, 0.71, 0.46, 0.39, 0.32, 0.57, 1.06, 1.38, 2.29, 2.99, 3.74, 4.39, 5.38, 7.25, 7.00, 4.57, 2.75, 1.48, 0.68, 0.39, 0.14 }; const char tract[] = "o";
constexpr int TractN = sizeof(A) / sizeof(A[0]);
const double l   = 0.397;
const double rho = 0.00114;
const double mu  = 0.000186;
const double c   = 35000;

constexpr double R(const int i) { return M_SQRTPI * std::sqrt(2. * rho * mu * om1) * std::pow(A[i], -1.5); }
constexpr double L(const int i) { return rho * l / (2. * A[i]); }
constexpr double C(const int i) { return l / (rho * c * c) * A[i]; }
constexpr double Ro(const int n) { return 128. * c * rho / (9. * M_PI * M_PI * A[n-1]); }
constexpr double Lo(const int n) { return  8. * rho / (3. * M_PI * std::sqrt(M_PI * A[n-1])); }

// y[] arrangement : [H1, V1, H2, V2, Ug, U1, ..., Un, P1, ..., Pn, Uo]
// int p[0] : current tract length
auto func = [](double t, const double y[], double dydt[], void *params) -> int
{
    const int n = TractN;
    const double h1 = Hmin > 0
                      ? (y[0] > 0 && y[0] < Hmin ? Hmin : y[0])
                      : y[0];
    const double v1 = y[1];
    const double h2 = Hmin > 0
                      ? (y[2] > 0 && y[2] < Hmin ? Hmin : y[2])
                      : y[2];
    const double v2 = y[3];
    const double Ug = y[4];
    const double *const U = y + 5;
    const double *const P = y + 5 + n;
    const double Uo = y[5 + n * 2];
    double *const dh1 = dydt + 0;
    double *const dv1 = dydt + 1;
    double *const dh2 = dydt + 2;
    double *const dv2 = dydt + 3;
    double *const dUg = dydt + 4;
    double *const dU = dydt + 5;
    double *const dP = dydt + 5 + n;
    double *const dUo = dydt + 5 + n * 2;

    const double x1 = h1 + xc;
    const double x2 = h2 + xc;
    const double s1 = ks1 * x1 * (1. + es1 * x1 * x1) + (h1 > 0 ? 0. : kh1 * h1 * (1. + eh1 * h1 * h1));
    const double s2 = ks2 * x2 * (1. + es2 * x2 * x2) + (h2 > 0 ? 0. : kh2 * h2 * (1. + eh2 * h2 * h2));
    const double r1 = 2. * m1 * om1 * (h1 > 0 ? xi1 : xi1 + 1.);
    const double r2 = 2. * m2 * om2 * (h2 > 0 ? xi2 : xi2 + 1.);

    const double Rbb = rho / (8. * lg * lg);
    const double Rv = 1.5 * mu / lg;
    const double Lg = rho / (2. * lg);
    const double Rc = 1.37 * Rbb * Ug / (h1 * h1);
    const double Rv1 = Rv * d1 / gsl_pow_3(h1);
    const double Lg1 = Lg * d1 / h1;
    const double R12 = Rbb * Ug * (1 / (h2 * h2) - 1 / (h1 * h1));
    const double Rv2 = Rv * d2 / gsl_pow_3(h2);
    const double Lg2 = Lg * d2 / h2;
    const double Re = -0.5 * Rbb * Ug / (h2 * h2);

    *dUg = (h1 <= 0 || h2 <= 0)
           ? -Ug / dt
           :  (Ps - P[0] - (Rc + Rv1 + R12 + Rv2 + Re) * Ug) / (Lg1 + Lg2);

    // need dUg
    const double dpi = -Rc * Ug;
    const double dp1 = -Rv1 * Ug - Lg1 * *dUg;
    const double dp12 = -R12 * Ug;
    const double dp2 = -Rv2 * Ug - Lg2 * *dUg;
//    const double dpo = Re * Ug;
    const double p1 = Ps + (h1 <= 0 || h2 <= 0 ? 0. : dpi + dp1 / 2);
    const double p2 = (h1 <= 0 ? 0. : Ps + (h2 <= 0 ? 0. : dpi + dp1 + dp12 + dp2 / 2));

    for (int i = 0; i < n; i++)
        dU[i] = (P[i] - P[i+1] - (R(i) + R(i+1)) * U[i]) / (L(i) + L(i+1));
    dU[n-1] = (P[n-1] - R(n-1) * U[n-1]) / L(n-1);
    *dUo = dU[n-1] - Ro(n) * Uo / Lo(n);

    dP[0] = (Ug - U[0]) / C(0);
    for (int i = 1; i < n; i++)
        dP[i] = (U[i-1] - U[i]) / C(i);

    *dh1 = v1;
    *dv1 = (p1 * lg * d1 - kc * (h1 - h2) - s1 - r1 * v1) / m1;
    *dh2 = v2;
    *dv2 = (p2 * lg * d2 - kc * (h2 - h1) - s2 - r2 * v2) / m2;

    return GSL_SUCCESS;
};

int main(int argc, char **argv)
{
    if (argc - 1 < 5)
    {
        std::cerr << "Too Few Arguments;" << std::endl;
        exit(EXIT_FAILURE);
    }

    const int d = 5 + TractN * 2 + 1;
    int p[1];
    gsl_odeiv_system sys = { func, NULL, d, p };

    gsl_odeiv_step *step = gsl_odeiv_step_alloc(StepType, d);

    double t = 0;
    double y[d] = {0};
    double dydt[d];
    double yerr[d];
    p[0] = sizeof(A) / sizeof(double);
    y[0] = std::atof(argv[1]);
    y[1] = std::atof(argv[2]);
    y[2] = std::atof(argv[3]);
    y[3] = std::atof(argv[4]);
    y[4] = std::atof(argv[5]);

    int m = (int)(T / dt);
    for (int i = 0; i <= m; i++)
    {
        t = dt * i;
        int status = gsl_odeiv_step_apply(step, t, dt, y, yerr, i == 0 ? NULL : dydt, dydt, &sys);
        if (status)
        {
            std::cerr << "failed : " << status << std::endl;
            exit(EXIT_FAILURE);
        }
        else if (gsl_isnan(yerr[0]))
        {
            std::cerr << "NaN error" << std::endl;
            exit(EXIT_FAILURE);
        }
        std::printf("%.6f %11.4e %11.4e %11.4e %11.4e %11.4e %11.4e %11.4e %11.4e %11.4e\n",
                    t,
                    y[0],
                    y[1],
                    y[2],
                    y[3],
                    2. * lg * y[0],
                    2. * lg * y[2],
                    y[4],
                    y[5 + TractN * 2],
                    dydt[5 + TractN * 2]
            );
    }

    return EXIT_SUCCESS;
}

https://twitter.com/suzukikojiro/status/535444806684520448