Metoda relaksacji (Richardsona)

Witam!

Mam do rozwiązania następujący układ równań metodą relaksacyjną.

Pierwsze z równań przekształciłem do następującej postaci:

Dla przypomnienia metoda relaksacyjna:

Problem polega na tym, że zbieżność tego algorytmu jest koszmarnie wolna (muszę wykonywać miliony iteracji) i nie wiem czy nie popełniłem gdzieś jakiegoś trudnego do zauważenia, a prostego do naprawienia błędu.

Poniżej zamieszczam mój kod:

#include <iostream>
#include <cmath>
#include <vector>
#include <iomanip>

constexpr unsigned N = 1001;
constexpr double h = 1.0 / (N - 1);
constexpr double h_squared = h * h;

void prepare_data(std::vector<double>& b) {
	b.resize(N, 0.0);
    b.front() = 1.0;
    b.back() = 1.0;
}

std::vector<double> richardson_iteration(std::vector<double>& b, double gamma, std::vector<double>& x) {
    for(unsigned k = 0; k < 9999999; ++k) {
        x[0] = x[0] + gamma * (b[0] - (1.0 * x[0]));
        for(unsigned i = 1; i < N - 1; ++i) {
            x[i] = x[i] + gamma * (b[i] - ((-1.0 * x[i - 1]) + ((2 + h_squared) * x[i]) + (-1.0 * x[i + 1])));
        }
        x[N - 1] = x[N - 1] + gamma * (b[N - 1] - (1.0 * x[N - 1]));
    }

	return x;
}

int main(int argc, char* argv[]) {
	std::vector<double> b;
	std::vector<double> q(N, 0.0);

	prepare_data(b);
    std::vector<double> x = richardson_iteration(b, 0.0000009999, q);

    for(auto e : x) {
        std::cout << std::fixed << std::setprecision(15) << std::endl;
    }

	return 0;
}

Pozdrawiam!

Nie Masz nigdzie ustawionego warunku przerwania przy jakiejś dokładności, ale zaglądając do środka algorytmu (print debug) widac, że szybko jest zbieżny do jedynki. Czy tak ma być? Na wikipedii pisze, że jak macierz spełnia pewne założenia to musi być zbieżny i to szybko (bo to to samo co GD).

No właśnie i tutaj jest kolejne pytanie. Generalnie dokładność jaką mam uzyskać to 10^(-10) natomiast nie wiem jak określić warunek stopu. Czy ma być to różnica norm tych dwóch wektorów mniejsza od jakiegoś epsilona czy raczej stopujemy kiedy maksymalna wartość różnic odpowiadających składowych jest mniejsza od epsilona?

Pierwsza i ostatnia składowa mają zmierzać do 1 i to się zgadza. Natomiast niepokoi mnie to co dzieje się dla np u_500 bo tam dzieje się... w zasadzie to nie wiem co.

Warunek stopu to to , że wektor x zmierza do jakiegoś punktu, czyli różnica kolejnych wektorków w iteracji dązy do zera, i U ciebie to sie szybko dzieje, środek do zera brzegi do jeden. W x[500] po prostu jest zero po jakimś czasie. Ustaw sobie, sprawdzanie co ileś iteracji czy róznica między dwoma wektorami jest mniejsza od dokładności (to 1000 działań, więc tzreba uważać)

Problem polega niestety na tym, że wiem jak ma wyglądać rozwiązanie i z pewnością środek nie dąży do zera. u_500 powinno wynosić około 0.886818892493490.

Ten algorytm lepiej się nie zbiegnie - wektor się za mało zmienia, około 10e-6 jest dokładność. Zrób to co zwykle w takich wypadkach, Sprawdź go dla jakis danych 2x2 i jak działa to takie jest życie, (może dla tych danych tak jest, mozna też metody optymalizcji, jak dla GD stosować) a jak nie to dupa debuging.