C++ Vec_DP类代码示例

void plotInterp(const Vec_DP &xdata, const Vec_DP &ydata,
        const Vec_DP &interpx, const Vec_DP &interpy) {
    int pointnum = 100;
    double start = xdata[0];
    double end = xdata[xdata.size() - 1];

    myplot::data_set container;
    myplot::plot_data realPoints("k+2");
    myplot::plot_data interpPoints("r.");
    myplot::plot_data curve("c-2");
    for (int i = 0; i < pointnum; i++) {
        double x = start + ((double) i * (end - start)) / (pointnum - 1);
        double y, dy;
        NR::polint(xdata, ydata, x, y, dy);
        curve.add_point(x, y);
    }
    for (int i = 0; i < xdata.size(); i++) {
        realPoints.add_point(xdata[i], ydata[i]);
    }
    for (int i = 0; i < interpx.size(); i++) {
        interpPoints.add_point(interpx[i], interpy[i]);
    }
    container.push_back(curve);
    container.push_back(realPoints);
    container.push_back(interpPoints);
    myplot::plot(container);
}

void generateNoisyData(Vec_DP& xdata, Vec_DP& ydata, Vec_DP& coeffs, double low, double high){
    for (int i = 0; i < coeffs.size(); i++) {
        coeffs[i] = rdm(low, high);
    }
    for (int i = 0; i < xdata.size(); i++) {
        ydata[i] = (coeffs[0]*pow(xdata[i],2) + coeffs[1]*xdata[i] + coeffs[2]) * rdm(0.8, 1.2);
    }
}

Vec_DP operator-(Vec_DP const& a, Vec_DP const& b) {
    if (a.size() != b.size())
        throw runtime_error("Dimensions mismatch in vector subtraction!");
    Vec_DP x(a.size());
    for (int i = 0; i < a.size(); i++) {
        x[i] = a[i] - b[i];
    }
    return x;
}

void spline_interp(Vec_DP const& xdata, Vec_DP const& ydata, double yp1, double ypn, Vec_DP const& x, Vec_DP& y) {
    Vec_DP y2(xdata.size());
    NR::spline(xdata, ydata, yp1, ypn, y2);
    for (int i = 0; i < x.size(); i++) {
        NR::splint(xdata, ydata, y2, x[i], y[i]);
    }

    double err, errx;
    get_error(x, y, err, errx);

    cout << "Maximum error in spline interpolation: " << err << ", at x = " << errx << endl;
}

void poly_interp(Vec_DP const& xdata, Vec_DP const& ydata, Vec_DP const& x, Vec_DP& y) {
    Vec_DP c(xdata.size());
    NR::polcoe(xdata, ydata, c);
    for (int i = 0; i < x.size(); i++) {
        y[i] = poly_val(x[i], c);
    }

    double err, errx;
    get_error(x, y, err, errx);

    cout << "Maximum error in polynomial interpolation: " << err << ", at x = " << errx << endl;
}

double interp_spline_periodic(const Vec_DP &xdata, const Vec_DP &ydata, double deriv,
        Vec_DP &interpx, Vec_DP &interpy, Vec_DP &error) {
    Vec_DP y2(xdata.size());
    NR::spline(xdata, ydata, deriv, deriv, y2);
    for (int i = 0; i < interpx.size(); i++) {
        double y, x = i * (M_PI / 200);
        NR::splint(xdata, ydata, y2, x, y);
        interpx[i] = x;
        interpy[i] = y;
        error[i] = abs(sin(x) - y);
    }
    return get_max_error(error);
}

void solve_system(Vec_DP xdata, Vec_DP ydata, Vec_DP zdata){
    Mat_DP A(xdata.size(), 3);
    Vec_DP b(xdata.size());
    for (int i = 0; i < xdata.size(); i++) {
        double x = xdata[i];
        double y = ydata[i];
        
        A[i][0] = x;
        A[i][1] = y;
        A[i][2] = 1;
        b[i] = zdata[i];
    }
    SVDsolve(A, b, coeffs);
}

void construct_data(Vec_DP &xdata, Vec_DP &ydata){
    for (int i = 0; i < xdata.size(); i++) {
        double x = i * (M_PI / 10);
        xdata[i] = x;
        ydata[i] = sin(x);
    }
}

double get_max_error(const Vec_DP &error){
    double max = 0;
    for (int i = 0; i < error.size(); i++) {
        max = (error[i] > max) ? error[i] : max;
    }
    return max;
}

double poly_val(double x, Vec_DP const& c) {
    double y = c[0];
    for (int n = 1; n < c.size(); n++) {
        y += c[n] * pow(x, n);
    }
    return y;
}

Vec_DP subst(Vec_DP const& w) {
    if (w.size() != 2)
        throw runtime_error("Wrong dimension for input in subst!");
    Vec_DP x(2);
    x[0] = -2 + 2.8 * pow(sin(w[0]), 2);
    x[1] = w[1];
    return x;
}

void constructData_B(Vec_DP &xdata, Vec_DP &ydata) {
    double ylow = -10;
    double yhigh = 10;
    for (int i = 0; i < xdata.size(); i++) {
        xdata[i] = i + 1;
        ydata[i] = rdm(ylow, yhigh);
    }
}

void interp(const Vec_DP &xdata, const Vec_DP &ydata,
        Vec_DP &interpx, Vec_DP &interpy, Vec_DP &error) {
    for (int i = 0; i < interpx.size(); i++) {
        double y, errory, x = 2 + i * 0.25;
        NR::polint(xdata, ydata, x, y, errory);
        interpx[i] = x;
        interpy[i] = y;
        error[i] = errory;
    }
}

void GSsolve_iter(Mat_DP &a, Vec_DP &b, Vec_DP &x, Vec_DP &xold) {
    int i, j, rows = a.nrows();
    if ((rows != a.ncols()) || (rows != b.size()) ||
            (rows != x.size()) || (rows != xold.size())) {
        cout << "Argument error in GaussSeidel_iter\n" << endl;
        abort();
    }

    double s;
    for (i = 0; i < rows; i++) {
        s = 0;
        for (j = 0; j < rows; j++) {
            s += (j <= i - 1) ? a[i][j] * x[j] :
                    (j >= i + 1) ? a[i][j] * xold[j] : 0;
        }
        xold[i] = x[i];
        x[i] = (b[i] - s) / a[i][i];
    }
}

double interp_polint(const Vec_DP &xdata, const Vec_DP &ydata,
        Vec_DP &interpx, Vec_DP &interpy, Vec_DP &error) {
    for (int i = 0; i < interpx.size(); i++) {
        double y, errory, x = i * (M_PI / 200);
        NR::polint(xdata, ydata, x, y, errory);
        interpx[i] = x;
        interpy[i] = y;
        error[i] = abs(sin(x) - y);
    }
    return get_max_error(error);
}