main.cpp - main - starlv/Fatigue_curve_estimation_cpp

main.cpp
143 строки5 KB
Agamirov L.V.
Add files via upload
22 май 2026, 12:27
22 май 2026, 12:27d2daaf4
#include "fatigue_curve.cpp"
#include <iostream>
#include <fstream>
#include <iomanip>
#include <vector>
#include <string>

using namespace std;

int main() {
    string ff = "fatigue";
    ifstream inp(ff + ".inp");
    ofstream out(ff + ".out");
    
    if (!inp.is_open()) {
        cerr << "Error: Cannot open input file " << ff << ".txt" << endl;
        return 1;
    }
    
    // Чтение данных из файла
    string st;
    int k, ku;
    double z;
    
    vector<double> cx;
    vector<int> ni;
    vector<double> lgN;
    vector<double> slgN;
    
    inp >> st;  // Levels_size
    inp >> k;
    inp >> st;  // Curve_type
    inp >> ku;
    inp >> st;  // Amplitudes
    for (int i = 0; i < k; i++) {
        inp >> z;
        cx.push_back(z);
    }
    inp >> st;  // Samples_size
    for (int i = 0; i < k; i++) {
        int kn;
        inp >> kn;
        ni.push_back(kn);
    }
    inp >> st;  // lgN_mean
    for (int i = 0; i < k; i++) {
        inp >> z;
        lgN.push_back(z);
    }
    inp >> st;  // SlgN
    for (int i = 0; i < k; i++) {
        inp >> z;
        slgN.push_back(z);
    }
    
    inp.close();
    
    // Вывод весов

    vector<double> w = computeWeights(ku, ni, lgN, slgN);

    out << "===== WEIGHTS =====" << endl;
    for (size_t i = 0; i < w.size(); ++i) {
        out << "Weight[" << i << "] = " << w[i] << endl;
    }
    
    // Оценка параметров кривой усталости
    double sigma_inf_guess=cx[k-1]/2.0;
    FatigueResult result = estimateFatigueCurve(ku, cx, ni, lgN, w, sigma_inf_guess);
    
    // Вывод результатов
    out << "\n========================================" << endl;
    out << "NONLINEAR LEAST SQUARES - Eigen LM" << endl;
    out << "========================================" << endl;
    if (ku == 0) out << "Model: lgN = (log10(C)-log10(σ-σ_inf))/m" << endl;
    if (ku == 1) out << "Model: N = (log10(C)-log10(σ-σ_inf))/m" << endl;
    
    out << "\nOptimization status: ";
    switch(result.optimization_status) {
        case 1: out << "Relative function tolerance achieved" << endl; break;
        case 2: out << "Relative parameter tolerance achieved" << endl; break;
        case 3: out << "Relative gradient tolerance achieved" << endl; break;
        case 4: out << "Max iterations reached" << endl; break;
        default: out << "Unknown (" << result.optimization_status << ")" << endl;
    }
    out << "Iterations: " << result.iterations << endl;
    
    out << "\nFinal parameters:" << endl;
    out << "sigma_inf = " << fixed << setprecision(2) << result.sigma_inf << " MPa" << endl;
    out << "C = " << scientific << setprecision(4) << result.C << endl;
    out << "m = " << fixed << setprecision(4) << result.m << endl;
    out << "Q (reduced chi-square) = " << fixed << setprecision(7) << result.Q << endl;
    
    // Ковариационная матрица
    out << "\nCovariance matrix:" << endl;
    out << "      [sigma_inf    log10(C)    m]" << endl;
    for (int i = 0; i < 3; ++i) {
        out << "Row" << i+1 << " ";
        for (int j = 0; j < 3; ++j) {
            out << scientific << setprecision(6) << setw(14) << result.covariance(i, j);
        }
        out << endl;
    }
    
    // Стандартные ошибки
    out << "\nStandard errors:" << endl;
    out << "se(sigma_inf) = " << scientific << sqrt(fabs(result.covariance(0, 0))) << endl;
    out << "se(log10(C))  = " << scientific << sqrt(fabs(result.covariance(1, 1))) << endl;
    out << "se(m)         = " << scientific << sqrt(fabs(result.covariance(2, 2))) << endl;
    
    out << "\n========================================" << endl;
    
    // Сравнение экспериментальных и расчетных значений
    out << "\nComparison:" << endl;
    if (ku == 0)
        out << "i\tσ(MPa)\t\tlgN_exp\t\tlgN_calc\tResidual" << endl;
    else
        out << "i\tσ(MPa)\t\tN_exp\t\tN_calc\t\tResidual" << endl;
    out << "------------------------------------------------------------------------------" << endl;
    
    for (size_t i = 0; i < cx.size(); ++i) {
        double pred, residual;
        if (ku == 0) {
            pred = pow(10, (log10(result.C) - log10(cx[i] - result.sigma_inf)) / result.m);
            residual = pred - lgN[i];
            out << i+1 << "\t" << fixed << setprecision(1) << cx[i] 
                << "\t\t" << setprecision(6) << lgN[i] 
                << "\t" << setprecision(6) << pred 
                << "\t" << scientific << setprecision(6) << residual << endl;
        } else { // ku == 1
            pred = (log10(result.C) - log10(cx[i] - result.sigma_inf)) / result.m;
            residual = pred - lgN[i];
            out << i+1 << "\t" << fixed << setprecision(1) << cx[i] 
                << "\t\t" << setprecision(6) << lgN[i] 
                << "\t" << setprecision(6) << pred 
                << "\t" << scientific << setprecision(6) << residual << endl;
        }
    }
    
    out.close();
    
    return 0;
}