Обратная задача 1 [Решена]

//
// Файл example1.txt
//
1 3

0    0 0    100  0 0
200  0 0    300  0 0
500  0 0    600  0 0
1000 0 0    1100 0 0

0.00530516
0.000265258
0.0000321525



//
// Файл Source.cpp
//

#include <iostream>
#include <cassert>
#include <fstream>
#include <vector>
#include <cmath>

using namespace std;

#define M_PI 3.1415926
#define EPS  0.0000001

struct Point3d {
    double x;
    double y;
    double z;
};

struct Range {
    Point3d start;
    Point3d end;
};

Range source; // Источник
vector<Range>  receivers;  // Приёмники
vector<double> potentials; // Разность потенциалов измеренное на каждом приёмнике

double calculatedI = 0.05;
vector<double> calculatedPotential;

constexpr double sigma = 0.1;
constexpr double goldenI = 1;




double distance(Point3d first, Point3d second) {
    return sqrt(pow(first.x - second.x, 2) +
                pow(first.y - second.y, 2) +
                pow(first.z - second.z, 2));
}

void calculateNextPotentials() {
    // Генерация синтетический данных
    for (size_t i = 0; i < receivers.size(); i++) {
        double U = calculatedI / (2 * M_PI * sigma) * ((1 / distance(source.end, receivers[i].start) - 1 / distance(source.start, receivers[i].start)) -
                                                       (1 / distance(source.end, receivers[i].end) - 1 / distance(source.start, receivers[i].end)));
        calculatedPotential[i] = U;
        std::cout << "calculatedPotential[" << i << "] = " << U << std::endl;
    }
}

// Минимизируемый функционал
double getF() {
    double F = 0;

    for (size_t i = 0; i < receivers.size(); i++) {
        F += 1 / potentials[i] * pow(potentials[i] - calculatedPotential[i], 2);
    }

    return F;
}

int main() {
    ifstream input("example1.txt");
    assert(input.is_open() && "Error while opening file");

    size_t sourcesCount;
    size_t receiversCount;

    input >> sourcesCount;
    input >> receiversCount;

    calculatedPotential.resize(receiversCount);

    // Считываем источник
    input >> source.start.x >> source.start.y >> source.start.z;
    input >> source.end.x   >> source.end.y   >> source.end.z;

    // Считываем приёмники
    for (size_t i = 0; i < receiversCount; i++) {
        Range receiver;

        input >> receiver.start.x >> receiver.start.y >> receiver.start.z;
        input >> receiver.end.x   >> receiver.end.y   >> receiver.end.z;

        receivers.emplace_back(receiver);
    }


    // Считываем разность потенциалов на каждом приёмнике
    for (size_t i = 0; i < receiversCount; i++) {
        double potential;

        input >> potential;

        potentials.emplace_back(potential);
    }
    vector<double> diff(receivers.size(), 0);
    calculateNextPotentials();


    // ---------------------- //
    while (getF() > EPS) {
        double elemA = 0; // Элемент СЛАУ
        double elemB = 0; // Правая часть уравнения
        for (size_t i = 0; i < receivers.size(); i++) {
            diff[i] = 1 / (2 * M_PI * sigma) * ((1 / distance(source.end, receivers[i].start) - 1 / distance(source.start, receivers[i].start)) -
                (1 / distance(source.end, receivers[i].end) - 1 / distance(source.start, receivers[i].end)));

            elemA += pow((1 / potentials[i]) * diff[i], 2);
            elemB -= pow(1 / potentials[i], 2) * diff[i] * (calculatedPotential[i] - potentials[i]);
        }
        double deltaI = elemB / elemA;
        calculatedI += deltaI;

        // Пересчитываем potenrial
        calculateNextPotentials();

        std::cout << calculatedI << std::endl;
    }
    system("pause");
    return 0;
}