-
Notifications
You must be signed in to change notification settings - Fork 0
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
Showing
17 changed files
with
162 additions
and
162 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -1,2 +1,2 @@ | ||
function next_step = AB2_next_step(f, D, x, y, A, B, h) | ||
next_step = y(1, :) + (h/2)*(3*D(2,:) - 1*D(1, :)); % predykcja | ||
next_step = y(1, :) + (h/2)*(3*D(2,:) - 1*D(1, :)); // prediction |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -1,7 +1,7 @@ | ||
function corr_step = M_corr_step(f, D, x, y, A, B, h, Corr_N) | ||
% x = x(i+1), y = y(i:i+1, :) | ||
corr_step = y(2, :); | ||
for c = 1:Corr_N % Zastosowanie korekcji Corr_N razy | ||
for c = 1:Corr_N % Applying correction Corr_N times | ||
DD = f(x, y(2, :), A, B); | ||
corr_step = y(1, :) + (h/2)*(DD + D(2, :)); | ||
end |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -1,2 +1,2 @@ | ||
function y = my_error(wiersz_gora,wiersz_dol) | ||
y = abs(wiersz_gora-wiersz_dol)/max(1,abs(wiersz_dol)); | ||
function y = my_error(upper_row, lower_row) | ||
y = abs(upper_row - lower_row) / max(1, abs(lower_row)); |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -1,29 +1,29 @@ | ||
function [Q, err, n_wierszy, wiersz_dol] = smartRM(g, a, b, tol, M0, Kmin, Kmax) | ||
% Autor: Łukasz Kryczka | ||
% Funkcja licząca całkę z funkcji g, na przedziale a,b, z tolerancją | ||
% błedu tol, w adaptacyjny sposób od M0 przedziałów z minimalnie | ||
% Kmin krokami ekstrapolacji i maksymalnie Kmax krokami ekstrapolacji. | ||
function [Q, err, num_rows, bottom_row] = smartRM(g, a, b, tol, M0, Kmin, Kmax) | ||
% Author: Łukasz Kryczka | ||
% Function to calculate the integral of function g, on the interval a,b, with error tolerance tol, | ||
% using an adaptive approach with M0 intervals, with a minimum of Kmin extrapolation steps | ||
% and a maximum of Kmax extrapolation steps. | ||
|
||
H = (b-a)/(M0+1); % liczymy zlozona calke trapezow dla M0+1 przedzialow | ||
wiersz_gora = H*(sum(g(a+H:H:b-H))+(g(a)+g(b))/2); | ||
H = (b-a)/(M0+1); % calculate composite trapezoidal integral for M0+1 intervals | ||
top_row = H*(sum(g(a+H:H:b-H))+(g(a)+g(b))/2); | ||
|
||
err = tol+1; % inicjalizujemy poczatkowa wartosc bledu | ||
n_wierszy=2; % zaczynamy ekstrapolacje od 2 wiersza | ||
while (err > tol || n_wierszy-1 <= Kmin) && n_wierszy-1 < Kmax+1 | ||
% liczymy pierwsze wartosci nowego wiersza T(i,0) | ||
wiersz_dol = zeros(1,n_wierszy); | ||
H = (b-a)/((M0+1)*2^(n_wierszy-1)); | ||
wiersz_dol(1) = H*(sum(g(a+H:H:b-H))+(g(a)+g(b))/2); | ||
err = tol+1; % initialize the initial error value | ||
num_rows = 2; % start extrapolation from the 2nd row | ||
while (err > tol || num_rows-1 <= Kmin) && num_rows-1 < Kmax+1 | ||
% calculate the first values of the new row T(i,0) | ||
bottom_row = zeros(1,num_rows); | ||
H = (b-a)/((M0+1)*2^(num_rows-1)); | ||
bottom_row(1) = H*(sum(g(a+H:H:b-H))+(g(a)+g(b))/2); | ||
|
||
% ekstrapolujemy kolejne wartosci T(i,k) | ||
for k = 2:n_wierszy | ||
wiersz_dol(k) = (4^(k-1)*wiersz_dol(k-1) ... | ||
-wiersz_gora(k-1))/(4^(k-1)-1); | ||
% extrapolate the next values T(i,k) | ||
for k = 2:num_rows | ||
bottom_row(k) = (4^(k-1)*bottom_row(k-1) ... | ||
-top_row(k-1))/(4^(k-1)-1); | ||
end | ||
% liczymy blad i aktualizujemy wartosci algorytmu | ||
err = my_error(wiersz_gora(end),wiersz_dol(end)); | ||
wiersz_gora = wiersz_dol; | ||
n_wierszy = 1 + n_wierszy; | ||
% calculate the error and update the algorithm values | ||
err = my_error(top_row(end),bottom_row(end)); | ||
top_row = bottom_row; | ||
num_rows = 1 + num_rows; | ||
end | ||
Q = wiersz_dol(end); % najdokladniejsza wartosc to ostatnia z wiersza | ||
n_wierszy = n_wierszy-2; % ilosc ekstrapolacji to n_wierszy-2 | ||
Q = bottom_row(end); % the most accurate value is the last one in the row | ||
num_rows = num_rows-2; % the number of extrapolations is num_rows-2 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -1,22 +1,22 @@ | ||
function wagi = baryc(wezly) | ||
% Autor: Łukasz Kryczka | ||
% Funkcja licząca wagi dla postaci barycentrycznej interpolacji | ||
% Lagrange'a | ||
function weights = baryc(nodes) | ||
% Author: Łukasz Kryczka | ||
% Function to calculate weights for barycentric Lagrange interpolation | ||
|
||
wagi = zeros(size(wezly)); | ||
if iscolumn(wezly) | ||
wezly=wezly'; | ||
weights = zeros(size(nodes)); | ||
if iscolumn(nodes) | ||
nodes = nodes'; | ||
end | ||
|
||
j = 0; | ||
for xj = wezly | ||
wspolczynnik = 1; | ||
for xi = wezly | ||
for xj = nodes | ||
coefficient = 1; | ||
for xi = nodes | ||
if (xj == xi) | ||
continue | ||
end | ||
wspolczynnik = wspolczynnik * (xj - xi); | ||
coefficient = coefficient * (xj - xi); | ||
end | ||
j = j + 1; | ||
wagi(j) = 1 ./ wspolczynnik; | ||
end | ||
weights(j) = 1 ./ coefficient; | ||
end | ||
} |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -1,19 +1,19 @@ | ||
function Y = ipbval(X, wezly, wagi, wartosci) | ||
% Autor: Łukasz Kryczka | ||
% Funkcja licząca wartość wielomianu interpolacyjnego Lagrang'a w | ||
% postaci barycentrycznej dla wektora X | ||
function Y = ipbval(X, nodes, weights, values) | ||
% Author: Łukasz Kryczka | ||
% Function to calculate the value of the Lagrange interpolation polynomial | ||
% in barycentric form for the vector X | ||
|
||
gj = 0; | ||
gj_f_xj = 0; | ||
|
||
for j = 1:length(wagi) | ||
gj = gj + wagi(j) ./ (X-wezly(j)); | ||
gj_f_xj = gj_f_xj + wagi(j) * wartosci(j) ./ (X-wezly(j)); | ||
for j = 1:length(weights) | ||
gj = gj + weights(j) ./ (X-nodes(j)); | ||
gj_f_xj = gj_f_xj + weights(j) * values(j) ./ (X-nodes(j)); | ||
end | ||
Y = gj_f_xj ./ gj; | ||
|
||
% x rowne wezla zamieniamy na ich oryginalne wartosci | ||
for i = find(ismember(X, wezly)) | ||
Y(i) = wartosci(wezly == X(i)); | ||
% Replace X equal to a node with its original value | ||
for i = find(ismember(X, nodes)) | ||
Y(i) = values(nodes == X(i)); | ||
end | ||
end |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -1,20 +1,20 @@ | ||
function wykwi(g,a,b,n) | ||
% Autor: Łukasz Kryczka | ||
% Funckja generująca wykresy dla wskaźnika na funkcję g, na przedziale | ||
% od a do b, dla n+1 węzłów. | ||
% Author: Łukasz Kryczka | ||
% Function generating plots for the function indicator g, on the interval | ||
% from a to b, for n+1 nodes. | ||
|
||
% wygenerowanie wezlow czebyszewa na przedziale od a do b | ||
% Generating Chebyshev nodes on the interval from a to b | ||
i = 0:n; | ||
wezly = (a + b) / 2 - ( (b-a)/2 ) .* cos( (2 .* i + 1)... | ||
nodes = (a + b) / 2 - ( (b-a)/2 ) .* cos( (2 .* i + 1)... | ||
.* pi ./ (2 * n + 2) ); | ||
X = linspace(a,b,2000000); | ||
|
||
% wykres funkcji | ||
% Plotting the function | ||
plot(X, g(X), color="b") | ||
hold on | ||
% wykres interpolacji funkcji | ||
plot(X, ipbval(X, wezly, baryc(wezly), g(wezly) ), color="r"); | ||
% wezly zaznaczone na zielono | ||
scatter(wezly, g(wezly), 'filled', 'green') | ||
% Plotting the interpolation of the function | ||
plot(X, ipbval(X, nodes, baryc(nodes), g(nodes) ), color="r"); | ||
% Nodes marked in green | ||
scatter(nodes, g(nodes), 'filled', 'green') | ||
hold off | ||
end |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -1,31 +1,31 @@ | ||
function composite_nodes = get_composite_nodes(a, b, n) | ||
% get_composite_nodes oblicza wektor przeskalowanych węzłów dla 3-punktowej | ||
% kwadratury Gaussa-Legendre'a dla złożonej kwadratury na przedziale [a, b] | ||
% get_composite_nodes calculates the vector of scaled nodes for a 3-point | ||
% Gauss-Legendre composite quadrature on the interval [a, b] | ||
% | ||
% Wejście: | ||
% a, b - granice przedziału całkowania, gdzie a < b, a i b rzeczywiste | ||
% n - liczba podprzedziałów, na które przedział [a, b] jest dzielony | ||
% liczba naturalna | ||
% Input: | ||
% a, b - integration limits, where a < b, a and b are real numbers | ||
% n - number of subintervals to divide the interval [a, b] into, | ||
% a natural number | ||
% | ||
% Wyjście: | ||
% composite_nodes - macierz węzłów kwadratury Gaussa-Legendre'a | ||
% Każda kolumna zawiera współczynniki dla odpowiedniego podprzedziału | ||
% Output: | ||
% composite_nodes - matrix of Gauss-Legendre quadrature nodes | ||
% Each column contains the coefficients for the corresponding subinterval | ||
% | ||
% Opis: | ||
% Funkcja get_composite_nodes działa przeskalowując i przesuwając | ||
% standardowe węzły kwadratury na każdy z n równych podprzedziałów | ||
% przedziału [a, b]. Wykorzystuje standardowe węzły 3-punktowej kwadratury | ||
% Gaussa-Legendre'a, skaluje je do rozmiaru każdego podprzedziału, | ||
% a następnie przesuwa do odpowiedniej lokalizacji. | ||
% Description: | ||
% The get_composite_nodes function scales and shifts the standard quadrature | ||
% nodes for each of the n equal subintervals of the interval [a, b]. | ||
% It uses the standard nodes of a 3-point Gauss-Legendre quadrature, | ||
% scales them to the size of each subinterval, | ||
% and then shifts them to the appropriate location. | ||
|
||
single_nodes=[-7.7459666924148337704e-01; | ||
0.0000000000000000000e+00; | ||
7.7459666924148337704e-01]; | ||
interval = (b-a)/(n); | ||
% Stworzenie wektora bazowego do skalowania i przesunięć | ||
% Create a base vector for scaling and shifting | ||
base = (a+interval/2):interval:b; | ||
% Powielenie single_nodes, aby pasowały do długości base | ||
% Replicate single_nodes to match the length of base | ||
nodes_replicated = repmat(single_nodes,1,length(base)); | ||
% Skalowanie i przesunięcie single_nodes | ||
% Scale and shift single_nodes | ||
composite_nodes = (interval/2)*nodes_replicated+base; | ||
end |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -1,38 +1,38 @@ | ||
function test_error() | ||
pause | ||
fprintf('Test dla wielomianu stopnia 6 - Start\n'); | ||
fprintf('Test for a polynomial of degree 6 - Start\n'); | ||
|
||
% Definiowanie przedziału i liczby podprzedziałów | ||
% Defining the interval and number of subintervals | ||
a = -1; b = 1; c = -1; d = 1; n = 1; m = 1; | ||
fprintf(['Testowanie dla przedziałów [%d, %d] i [%d, %d] z liczbą ' ... | ||
'podprzedziałów n = %d, m = %d.\n'], a, b, c, d, n, m); | ||
fprintf(['Testing for intervals [%d, %d] and [%d, %d] with number ' ... | ||
'of subintervals n = %d, m = %d.\n'], a, b, c, d, n, m); | ||
|
||
% Funkcja szóstego stopnia | ||
% Sixth degree function | ||
f = @(x, y) 1.2*x.^6 - y.^6 + 0.8*x.^5.*y; | ||
fprintf('f(x, y) = 1.2*x.^6 - y.^6 + 0.8*x.^5.*y\n'); | ||
|
||
% Oczekiwana wartość całki | ||
% Expected integral value | ||
expected_value = 4/35; | ||
fprintf('Oczekiwana wartość całki dla funkcji rzędu 6: %.2f\n', expected_value); | ||
fprintf('Metoda oparta na 3 węzłach, więc błąd powinien być substancjalny\n') | ||
fprintf('Expected integral value for a degree 6 function: %.2f\n', expected_value); | ||
fprintf('Method based on 3 nodes, so the error should be substantial\n') | ||
|
||
% (d-c)d1 + (b-a)d2 (c1 + d2) / m^6 | ||
% (d-c)d1 + (b-a)d2 (c1 + d2) / 2^6 | ||
% d1 = c1/m^2n | ||
% d2 = d1/n^2n | ||
% Przybliżona wartość całki obliczona przez funkcję | ||
approx_value = P1Z30_LKR_CDIGL(f, a, b, c, d, n, m); | ||
fprintf('Przybliżona wartość całki: %.2f\n', approx_value); | ||
% Błąd bezwzględny | ||
% Approximate integral value calculated by the function | ||
approx_value = double_integral_gauss_legendre(f, a, b, c, d, n, m); | ||
fprintf('Approximate integral value: %.2f\n', approx_value); | ||
% Absolute error | ||
error = abs(approx_value - expected_value); | ||
error_old = error; | ||
fprintf('Błąd bezwzględny: %.2e\n', error); | ||
fprintf('Absolute error: %.2e\n', error); | ||
n = 2; m = 2; | ||
approx_value = P1Z30_LKR_CDIGL(f, a, b, c, d, n, m); | ||
fprintf('Przybliżona wartość całki: %.2f\n', approx_value); | ||
fprintf('Approximate integral value: %.2f\n', approx_value); | ||
|
||
% Błąd bezwzględny | ||
% Absolute error | ||
error = abs(approx_value - expected_value); | ||
fprintf('Błąd bezwzględny: %.2e\n', error); | ||
fprintf('Absolute error: %.2e\n', error); | ||
|
||
fprintf('is it equal %d %d\n', error, error_old/(2^6)) |
24 changes: 12 additions & 12 deletions
24
Gauss-Legendre-2D-Integration/tests_composite_nodes/test_composite_nodes_basic.m
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -1,29 +1,29 @@ | ||
function test_composite_nodes_basic() | ||
pause | ||
disp('Test: Podstawowa funkcjonalność - Start'); | ||
disp('Test: Basic functionality - Start'); | ||
|
||
% Definiowanie przedziału i liczby podprzedziałów | ||
% Defining the interval and number of subintervals | ||
a = 0; b = 3; n = 3; | ||
disp(['Testowanie dla przedziału [', num2str(a), ', ', num2str(b), ... | ||
'] z liczbą podprzedziałów n = ', num2str(n)]); | ||
disp(['Testing for the interval [', num2str(a), ', ', num2str(b), ... | ||
'] with the number of subintervals n = ', num2str(n)]); | ||
|
||
% Oczekiwany wynik obliczony ręcznie | ||
% Expected result calculated manually | ||
expected_output = [1 1]; | ||
fprintf('Oczekiwany wynik:\n'); | ||
fprintf('Expected result:\n'); | ||
disp(expected_output); | ||
fprintf('\n'); | ||
|
||
% Pobieranie rzeczywistego wyniku z funkcji | ||
% Getting the actual result from the function | ||
actual_output = get_composite_nodes(a, b, n); | ||
fprintf('Rzeczywisty wynik:\n'); | ||
fprintf('Actual result:\n'); | ||
disp(actual_output); | ||
fprintf('\n'); | ||
|
||
% Sprawdzenie, czy rzeczywisty wynik zgadza się z oczekiwanym | ||
% Tolerancja jest ustawiona, aby uwzględnić błędy arytmetyki zmiennoprzecinkowej | ||
% Checking if the actual result matches the expected result | ||
% Tolerance is set to account for floating-point arithmetic errors | ||
difference = abs(actual_output - expected_output); | ||
fprintf('Różnica: ['); | ||
fprintf('Difference: ['); | ||
fprintf(' %.2e ', difference); | ||
fprintf(']\n'); | ||
fprintf('Test: Podstawowa funkcjonalność - Koniec\n'); | ||
fprintf('Test: Basic functionality - End\n'); | ||
end |
Oops, something went wrong.