Санкт-Петербурзький університет шляхів сполучення
Кафедра «Прикладна математика»
ЗВІТ
По виконанню курсових робіт
Предмет «Чисельні методи»
«Застосування чисельних методів для вирішення Рівнянь з приватними похідними»
Санкт-Петербург 2008р.
Лабораторна робота N1 "Інтерполяція алгебраїчними многочленами"
Для вирішення завдання локального інтерполяції алгебраїчними многочленами в системі MATLAB призначені функції polyfit (POLYnomial FITting - апроксимація многочленом) і polyval (POLYnomial VALue - значення многочлена).
Функція polyfit (X, Y, n) знаходить коефіцієнти многочлена ступеня n, побудованого за даними вектора Х, який апроксимує дані вектора Y в сенсі найменшого квадрата відхилення. Якщо число елементів векторів X і Y одно n +1, то функція polyfit (X, Y, n) вирішує задачу інтерполяції многочленом ступеню n.
Функція polyval (P, z) обчислює значення полінома, коефіцієнти якого є елементами вектора P, від аргументу z. Якщо z - вектор або матриця, то поліном обчислюється у всіх точках z.
Скористаємося зазначеними функціями системи MATLAB для вирішення задачі локального інтерполяції алгебраїчними многочленами функції, заданої таблицею своїх значень
X
0.0
1.0
2.0
3.0
4.0
Y
1.0
1.8
2.2
1.4
1.0
і обчислення її наближеного значення в точці x * = 2.2.
Задача 1 (завдання локального інтерполяції многочленами)
Побудувати інтерполяційні многочлени 1-ої, 2-ий і 3-го ступеня.
Обчислити їх значення при x = x *.
Записати багаточлени в канонічній формі і побудувати їх графіки.
Рішення завдання засобами системи MATLAB:
X = [0.0000 0.5000 1.0000 1.5000 2.0000 2.5000];
Y = [0.0378 0.0653 0.3789 1.0353 0.5172 0.9765];
xzv = 1.61;
P1 = polyfit (X (4:5), Y (4:5), 1) Коефіцієнти многочлена P1
P2 = polyfit (X (3:5), Y (3:5), 2) Коефіцієнти многочлена P2
P3 = polyfit (X (3:6), Y (3:6), 3) Коефіцієнти многочлена P3
Отримані таким чином коефіцієнти інтерполяційних многочленів і значення цих многочленів при x = x *:
P1 = -1.0362 2.5896
P2 = -2.3490 7.1853 -4.4574
P3 = 2.8692 -15.2604 25.8351 -13.0650
z1 = 0.9213
z2 = 1.0221
z3 = 0.9470
многочлени P1, P2, P3
P 1 = -1.0362 * X + 2.5896
P 2 = -2.3490 * X 2 + 7.1853 * X + -4.4574
P 3 = 2.8692 * X 3 -15.2604 * X 2 + 25.8351 + -13.0650
Для побудови графіків інтерполяційних многочленів слід створити вектори xi1, xi2, xi3, що моделюють інтервали (X (3): X (4)), (X (2): X (4)), (X (2): X (5) ), відповідно, і обчислити значення многочленів P1, P2, P3 для елементів векторів xi1, xi2, xi3, відповідно:
xi1 = X (4): 0.05: X (5);
xi2 = X (3): 0.05: X (5);
xi3 = X (3): 0.05: X (6);
y1 = polyval (P1, xi1);
y2 = polyval (P2, xi2);
y3 = polyval (P3, xi3);
plot (X, Y, '* k', xi1, y1, xi2, y2, xi3, y3); grid
Інтерполяція нелінійною функцією Y = A * exp (- B * X)
y_l = log (Y)
Pu = polyfit (X (4:5), y_l (4:5), 1)
z_l = (exp (Pu (2)) * exp (Pu (1) * xzv))
Y = 8.3040 * exp (-1.3880 * X)
Функція plot із зазначеними аргументами будує табличні значення функції чорними зірочками ('* k'), а також графіки многочленів P1 (по векторах xi1 і y1), P2 (по векторах xi2 і y2) і P3 (по векторах xi3 і y3), і функцією Y = A * exp (- B * X), відповідно синьою, червоною і зеленою кривими.
plot (X, Y, '* k', xi1, y1, xi2, y2, xi3, y3, xi1, exp (Pu (2)) * exp (Pu (1) * xi1)); grid
Оцінка похибки інтерполяції
При оцінці похибки розв'язку задачі інтерполяції у точці x * за похибка epsk інтерполяційного многочлена ступеня k приймається модуль різниці значень цього многочлена і многочлена ступеня k +1 в точці x *.
За допомогою вже отриманих значень ми можемо оцінити похибки інтерполяційних многочленів P1 і P2 в точці x *, використовуючи функцію abs системи MATLAB для обчислення модуля:
eps1 = abs (z1-z2)
eps1 = 0.1008
eps2 = abs (z2-z3)
eps2 = 0.0751
Для оцінки похибки многочлена P3 необхідно попередньо обчислити
значення z4 = P4 (x *), а потім - eps3.
P4 = polyfit (X, Y, 4); z4 = polyval (P4, xzv);
eps3 = abs (z4-z3)
eps3 = 0.1450
«Побудова сплайна»
X = [0.0000 0.5000 1.0000 1.5000 2.0000 2.5000];
Y = [0.0378 0.0653 0.3789 1.0353 0.5172 0.9765];
cs = spline (X, [0 Y 0]);
xx = linspace (0,2.5);
plot (X, Y, '* m', xx, ppval (cs, xx), '-k');
h = 0.5
esstestvennii spline
A = [4 2 0 0 0 0
1 4 1 0 0 0
0 1 4 1 0 0
0 0 1 4 1 0
0 0 0 1 4 1
0 0 0 0 2 4]
B = [6 * (Y (2)-Y (1)) / h 0 0 0 0 6 * (Y (length (Y))-Y (length (Y) -1)) / h]
for i = 2: (length (Y) -1)
B (i) = (3 / h) * (Y (i +1)-Y (i-1))
end
S = inv (A) * B '
otsutstvie uzla
A1 = [1 0 -1 0 0 0
1 4 1 0 0 0
0 1 4 1 0 0
0 0 1 4 1 0
0 0 0 1 4 1
0 0 0 1 0 -1]
B1 = [2 * (2 * Y (2)-Y (1)-Y (3)) / h 0 0 0 0 2 * (2 * Y (length (Y) -1)-Y (length (Y) )-Y (length (Y) -2)) / h]
for i = 2: (length (Y) -1)
B1 (i) = (3 / h) * (Y (i +1)-Y (i-1))
end
S1 = inv (A1) * B1 '
c1 = spline (X, [S (2) YS (5)]);
x1 = linspace (0,2.5,101);
c2 = spline (X, [S1 (2) Y S1 (5)]);
x2 = linspace (0,2.5,101);
plot (X, Y, 'ob', xx, ppval (cs, xx ),'-', x1, ppval (c1, x1 ),'*', x2, ppval (c2, x2 ),'^', xx , spline (X, Y, xx));
h = 0.5000
A =
4 2 0 0 0 0
1 4 1 0 0 0
0 1 4 1 0 0
0 0 1 4 1 0
0 0 0 1 4 1
0 0 0 0 2 4
B = 0.3300 0 0 0 0 5.5116
B = 0.3300 2.0466 0 0 0 5.5116
B = 0.3300 2.0466 5.8200 0 0 5.5116
B = 0.3300 2.0466 5.8200 0.8298 0 5.5116
B = 0.3300 2.0466 5.8200 0.8298 -0.3528 5.5116
S =
0.0052
0.1546
1.4230
-0.0266
-0.4869
1.6213
A1 =
1 0 -1 0 0 0
1 4 1 0 0 0
0 1 4 1 0 0
0 0 1 4 1 0
0 0 0 1 4 1
0 0 0 1 0 -1
B1 = -1.1444 0 0 0 0 -3.9096
B1 = -1.1444 2.0466 0 0 0 -3.9096
B1 = -1.1444 2.0466 5.8200 0 0 -3.9096
B1 = -1.1444 2.0466 5.8200 0.8298 0 -3.9096
B1 = -1.1444 2.0466 5.8200 0.8298 -0.3528 -3.9096
S1 =
0.2496
0.1008
1.3940
0.1433
-1.1372
4.0529
Лабораторна робота N2 "Побудова алгебраїчних многочленів найкращого середньоквадратичного наближення"
X = [0.0000 0.5000 1.0000 1.5000 2.0000 2.5000];
Y = [0.0378 0.0653 0.3789 1.0353 0.5172 0.9765];
n = length (X)
TABL = [X, sum (X); Y, sum (Y );...
X. ^ 2, sum (X. ^ 2 );...
X. * Y, sum (X. * Y );...
X. * X. * Y, sum (X. * X. * Y );...
X. ^ 3, sum (X. ^ 3); X. ^ 4, sum (X. ^ 4)];
TABL = TABL '
XYX ^ 2 X * YX ^ 2 * YX ^ 3 X ^ 4
0 0.0378 0 0 0 0 0
0.5000 0.0653 0.2500 0.0326 0.0163 0.1250 0.0625
1.0000 0.3789 1.0000 0.3789 0.3789 1.0000 1.0000
1.5000 1.0353 2.2500 1.5530 2.3294 3.3750 5.0625
2.0000 0.5172 4.0000 1.0344 2.0688 8.0000 16.0000
2.5000 0.9765 6.2500 2.4413 6.1031 15.6250 39.0625
7.5000 3.0110 13.7500 5.4402 10.8966 28.1250 61.1875 - Сума
За даними таблиці запишемо і вирішимо нормальну систему МНК-методу:
1) дл многочлена першого ступеня
S1 = [n, TABL (7,1); TABL (7,1) TABL (7,3)] матриця коефіцієнтів
T1 = [TABL (7,2); TABL (7,4)] вектор правих частин
coef1 = S1 \ T1 рішення нормальної системи МНК
A1 = coef1 (2); B1 = coef1 (1); коефіцієнти многочлена 1-го ступеня
S1 =
6.0000 7.5000
7.5000 13.7500
T1 =
3.0110
5.4402
coef1 =
0.0229
0.3832
2) дл многочлена другого ступеня
S2 = [n TABL (7,1) TABL (7,3); TABL (7,1) TABL (7,3) TABL (7,6); TABL (7,3) TABL (7,6) TABL ( 7,7)] матриця коефіцієнтів
T2 = [TABL (7,2); TABL (7,4); TABL (7,5)] вектор правих частин
coef2 = S2 \ T2 рішення нормальної системи МНК
A2 = coef2 (3); B2 = coef2 (2); C2 = coef2 (1); коефіцієнти многочлена 2-го ступеня
S2 =
6.0000 7.5000 13.7500
7.5000 13.7500 28.1250
13.7500 28.1250 61.1875
T2 =
3.0110
5.4402
10.8966
coef2 =
-0.0466
0.5917
-0.0834
Для побудови графіків функцій y1 = A1 * x + B1 і y2 = A2 * x ^ 2 + B2 * x + C2 зі знайденими коефіцієнтами задамо допоміжний вектор абсциси xi, а потім обчислимо елементи векторів g1 = A1 * xi + B1 і g2 = A2 * xi ^ 2 + B2 * xi + C2:
h = 0.05;
xi = min (X): h: max (X);
g1 = A1 * xi + B1;
g2 = A2 * xi. ^ 2 + B2 * xi + C2;
plot (X, Y, '* k', xi, g1, xi, g2); grid
coef1 = polyfit (X, Y, 1) коефіцієнти многочлена першого ступеня
coef2 = polyfit (X, Y, 2) коефіцієнти многочлена другого ступеня
coef1 = 0.3832 0.0229
coef2 = -0.0834 0.5917 -0.0466
Для побудови графіків задамо допоміжний вектор абсциси xi, а потім c допомогою функції polyval обчислимо елементи векторів g1 і g2:
xi = min (X): 0.1: max (X);
g1 = polyval (coef1, xi);
g2 = polyval (coef2, xi);
plot (X, Y, '* k', xi, g1, xi, g2); grid
Очевидно, що побудовані таким способом графіки співпадуть з отриманими раніше.
Для того, щоб визначити величину середньоквадратичного відхилення, обчислимо суми квадратів відхилень g1 (x) і g2 (x) від таблично заданої функції у вузлах таблиці X а потім
G1 = polyval (coef1, X);
G2 = polyval (coef2, X);
delt1 = sum ((Y-G1). ^ 2); delt1 = sqrt (delt1 / 5)
delt2 = sum ((Y-G2). ^ 2); delt2 = sqrt (delt2 / 5)
Останні два рядки можна замінити іншими, якщо використовувати функцію mean, вичіслющую середнє значення:
delt1 = mean (sum ((Y-G1). ^ 2))
delt2 = mean (sum ((Y-G2). ^ 2))
delt1 = 0.2403
delt2 = 0.2335
delt1 = 0.2888
delt2 = 0.2725
Для нелінійної
X = [0.0000 0.5000 1.0000 1.5000 2.0000 2.5000];
Y = [0.0378 0.0653 0.3789 1.0353 0.5172 0.9765]
Y_o = Y
Y = 1. / (Exp (Y))
n = length (X)
TABL = [X, sum (X); Y, sum (Y );... означає перенесення рядка
X. ^ 2, sum (X. ^ 2 );...
X. * Y, sum (X. * Y );...
X. * X. * Y, sum (X. * X. * Y );...
X. ^ 3, sum (X. ^ 3); X. ^ 4, sum (X. ^ 4)];
TABL = TABL '
За даними таблиці запишемо і вирішимо нормальну систему МНК-методу:
2) дл многочлена другого ступеня
S2 = [n TABL (7,1) TABL (7,3); TABL (7,1) TABL (7,3) TABL (7,6); TABL (7,3) TABL (7,6) TABL ( 7,7)] матриця коефіцієнтів
T2 = [TABL (7,2); TABL (7,4); TABL (7,5)] вектор правих частин coef2 = S2 \ T2 рішення нормальної системи МНК
A2 = coef2 (3); B2 = coef2 (2); C2 = coef2 (1); коефіцієнти многочлена 2-го ступеня
Дл побудови графіків функції y2 = A2 * x ^ 2 + B2 * x + C2 зі знайденими коефіцієнтами задамо допоміжний вектор абсциси xi, а потім обчислимо елементи векторів g1 = A1 * xi + B1 і g2 = A2 * xi ^ 2 + B2 * xi + C2:
h = 0.05;
xi = min (X): h: max (X);
g2 = log (1. / (A2 * xi. ^ 2 + B2 * xi + C2));
plot (X, Y_o, '* k', xi, g2); grid
coef2 = polyfit (X, Y, 2) коефіцієнти многочлена другого ступеня
Для побудови графіків задамо допоміжний вектор абсциси xi, а потім c допомогою функції polyval обчислимо елементи векторів g1 і g2:
pause;
xi = min (X): 0.1: max (X);
g2 = polyval (coef2, xi);
plot (X, Y_o, '* k', xi, log (1./g2)); grid
Очевидно, що побудовані таким способом графіки співпадуть з отриманими раніше.
Для того, щоб визначити величину середньоквадратичного відхилення, обчислимо суми квадратів відхилень g1 (x) і g2 (x) від таблично заданої функції у вузлах таблиці X а потім
G2 = polyval (coef2, X);
delt2 = sum ((Y-G2). ^ 2); delt2 = sqrt (delt2 / 5)
Останні два рядки можна замінити іншими, якщо використовувати функцію mean, яка обчислює середнє значення:
delt2 = mean (sum ((Y-G2). ^ 2))
Y = 0.0378 0.0653 0.3789 1.0353 0.5172 0.9765
Y_o = 0.0378 0.0653 0.3789 1.0353 0.5172 0.9765
Y = 0.9629 0.9368 0.6846 0.3551 0.5962 0.3766
n = 6
TABL =
0 0.9629 0 0 0 0 0
0.5000 0.9368 0.2500 0.4684 0.2342 0.1250 0.0625
1.0000 0.6846 1.0000 0.6846 0.6846 1.0000 1.0000
1.5000 0.3551 2.2500 0.5327 0.7990 3.3750 5.0625
2.0000 0.5962 4.0000 1.1924 2.3848 8.0000 16.0000
2.5000 0.3766 6.2500 0.9416 2.3539 15.6250 39.0625
7.5000 3.9122 13.7500 3.8196 6.4565 28.1250 61.1875
S2 =
6.0000 7.5000 13.7500
7.5000 13.7500 28.1250
13.7500 28.1250 61.1875
T2 =
3.9122
3.8196
6.4565
coef2 =
1.0178
-0.4243
0.0718
coef2 =
0.0718 -0.4243 1.0178
delt2 =
0.1199
delt2 =
0.0719
Чисельні методи розв'язання задачі Коші для звичайних диференціальних рівнянь
Ейлера явна
function dy = func (x, y)
dy = 2 * x * y
clear
X = [0.00000 0.10000 0.20000 0.30000 0.40000 0.50000];
Y = exp ((X). ^ 2);
Y0 = input ('Значення функції в точці 0 =');
Y_n1 = Y0;
for n = 1: length (X) -1
Y_n1 = Y_n1 +0.1 * 2 * X (n) * Y_n1;
Y_n (n) = Y_n1;
end
X1 = 0.00000:0.01:0.50000;
Y_sot = Y0;
for n = 1: length (X1)
Y_sot = Y_sot +0.01 * 2 * X1 (n) * Y_sot;
Y_sto (n) = Y_sot;
end
X2 = 0.00000:0.001:0.50000;
Y_tys = Y0;
for n = 1: length (X2)
Y_tys = Y_tys +0.001 * 2 * X2 (n) * Y_tys;
Y_ts (n) = Y_tys;
end
disp ('XY h = 0.1 h = 0.01 h = 0.001')
G1 = Y_sto (10:10: end);
TABL = [X, Y; Y0, Y_n; ...
Y_sto (1), G1; ...
Y_ts (1), Y_ts (100:100: end)];
TABL = TABL '; disp (TABL)
Значення функції в точці 0 = 1
X Y h = 0.1 h = 0.01 h = 0.001
0 1.0000 1.0000 1.0000 1.0000
0.1000 1.0101 1.0000 1.0090 1.0099
0.2000 1.0408 1.0200 1.0387 1.0406
0.3000 1.0942 1.0608 1.0907 1.0938
0.4000 1.1735 1.1244 1.1683 1.1730
0.5000 1.2840 1.2144 1.2766 1.2833
Симетрична
clear
X = [0.00000 0.10000 0.20000 0.30000 0.40000 0.50000];
Y = exp ((X). ^ 2);
Y0 = input ('Значення функції в точці 0 =');
Y_n1 = Y0;
for n = 1: length (X) -1
Y_n1 = Y_n1 * (1 +0.1 * X (n)) / (1-0.1 * (X (n) +0.1));
Y_n (n) = Y_n1;
end
X1 = 0.00000:0.01:0.50000;
Y_sot = Y0;
for n = 1: length (X1) -1
Y_sot = Y_sot * (1 +0.01 * X1 (n)) / (1-0.01 * (X1 (n) +0.01));
Y_sto (n) = Y_sot;
end
X2 = 0.00000:0.001:0.50000;
Y_tys = Y0;
for n = 1: length (X2) -1
Y_tys = Y_tys * (1 +0.001 * X2 (n)) / (1-0.001 * (X2 (n) +0.001));
Y_ts (n) = Y_tys;
end
disp ('XY h = 0.1 h = 0.01 h = 0.001')
G1 = Y_sto (10:10: end);
TABL = [X, Y; Y0, Y_n; ...
Y_sto (1), G1; ...
Y_ts (1), Y_ts (100:100: end)];
TABL = TABL '; disp (TABL)
Значення функції в точці 0 = 1
X Y h = 0.1 h = 0.01 h = 0.001
0 1.0000 1.0000 1.0001 1.0000
0.1000 1.0101 1.0101 1.0101 1.0101
0.2000 1.0408 1.0410 1.0408 1.0408
0.3000 1.0942 1.0947 1.0942 1.0942
0.4000 1.1735 1.1745 1.1735 1.1735
0.5000 1.2840 1.2858 1.2840 1.2840
Ейлера неявна
clc
clear all
h_1 = 0.1;
X = 0: h_1: 0.5;
Y = exp (X. ^ 2);
Yn = Y (1);
Y2 = Yn + h_1 * 2 * X (1);
або Y2 = input ('Введіть значення Yn в точці X = 0 =')
y_1 (1) = Y2;
for i = 1: (length (Y) -1)
y_1 (i +1) = y_1 (i) / (1-h_1 * 2 * X (i +1));
end
h_2 = 0.01;
X_2 = 0: h_2: 0.5;
Y_2 = exp (X_2. ^ 2);
Y2 = Yn + h_2 * 2 * X (1);
y_2 (1) = Y2;
for i = 1: (length (Y_2) -1)
y_2 (i +1) = y_2 (i) / (1-h_2 * 2 * X_2 (i +1));
end
h_3 = 0.001;
X_3 = 0: h_3: 0.5;
Y_3 = exp (X_3. ^ 2);
Y2 = Yn + h_3 * 2 * X (1);
y_3 (1) = Y2;
for i = 1: (length (Y_3) -1)
y_3 (i +1) = y_3 (i) / (1-h_3 * 2 * X_3 (i +1));
end
for k = 1:5
r1 (k) = y_2 (k * 10 +1);
r2 (k) = y_3 (k * 100 +1);
end
TABL = [X, Y; ... ... означає перенесення рядка
y_1; ...
y_2 (1), r1; ...
y_3 (1), r2];
TABL = TABL '
plot (X, Y ,'-', X, y_1, X, [y_2 (1), r1], X, [y_3 (1), r2])
f = ode23 ('y_1', [0 0.5], 1)
TABL =
0 1.0000 1.0000 1.0000 1.0000
0.1000 1.0101 1.0204 1.0111 1.0102
0.2000 1.0408 1.0629 1.0430 1.0410
0.3000 1.0942 1.1308 1.0977 1.0945
0.4000 1.1735 1.2291 1.1787 1.1740
0.5000 1.2840 1.3657 1.2916 1.2848
Завдання Коші
function [output_args] = koshi (input_args)
KOSHI Summary of this function goes here
Detailed explanation goes here
tspan = [0,1];
y0 = [1.421,1];
[T, y] = ode45 (@ F, tspan, y0);
ode45 (@ F, tspan, y0);
hold on
plot ([0 1], [1 1])
Підбір Альфа методом січних
a = 1;
y0 = [1, a];
tspan = [0,1];
psi_old = a-1;
a_old = 0.5;
i = 1;
eps = 1;
while (eps> = 0.000001) & (i <10000)
[T, y] = ode45 (@ F, tspan, y0);
ode45 (@ F, tspan, y0)
psi = y (end, 2) -1;
a_new = a-psi * (a-a_old) / (psi-psi_old)
eps = abs (psi);
a_old = a;
a = a_new;
y0 = [1, a];
psi_old = psi
i = i + 1;
end
i
a_new = 0.5000
psi_old = -0.3935
a_new = 1.4655
psi_old = -0.8161
a_new = 1.4184
psi_old = 0.0419
a_new = 1.4215
psi_old = -0.0030
a_new = 1.4215
psi_old =-4.1359e-006
a_new = 1.4215
psi_old = 4.2046e-010
i = 7
Генерація матриці 10 * 10 з елементів рівномірно розподілених 1 .. 10
function [output_args] = ravnomern10_10_1_10 (input_args)
RAVNOMERN10_10_1_10 Summary of this function goes here
Detailed explanation goes here
round (rand (10,10) * 9 +1)
8 2 7 7 5 3 8 9 4 2
9 10 1 1 4 7 3 3 8 1
2 10 9 3 8 7 6 8 6 6
9 5 9 1 8 2 7 3 6 8
7 8 7 2 3 2 9 9 9 9
2 2 8 8 5 5 10 4 4 2
4 5 8 7 5 10 6 3 8 6
6 9 5 4 7 4 2 3 8 5
10 8 7 10 7 6 2 7 4 1
10 10 3 1 8 3 3 5 6 4
Рішення краєвої задачі методом сіток для УЧП.
e = 0.01;
h = sqrt (e);
x = 0: h: 1;
y = 0: h: 1;
v = ones (11,11);
v (1,:) = 0;
v (end,:) = 1;
v (:, 1) = (0: h: 1) ';
v (:, end) = (0: h: 1) ';
v = v. * ((1 * 9 + sum (0: h: 1) + sum (0: h: 1)) / 40)
v (1,:) = 0;
v (end,:) = 1;
v (:, 1) = (0: h: 1) ';
v (:, end) = (0: h: 1) ';
surf (v);
d = e +1;
i = 1
while d> e & i <100
v1 = v;
v1 (1:10,:) = v1 (2:11,:);
v1 (11,:) = v (1,:);
v2 = v;
v2 (2:11,:) = v2 (1:10,:);
v2 (1,:) = v (11,:);
v3 = v;
v3 (:, 1:10) = v3 (:, 2:11);
v3 (:, 11) = v (:, 1);
v4 = v;
v4 (:, 2:11) = v4 (:, 1:10);
v4 (:, 1) = v (:, 11);
v_new = (v1 + v2 + v3 + v4) / 4;
d = max (max (abs (v (2: end-1, 2: end-1)-v_new (2: end-1, 2: end-1))))
v = v_new;
v (1,:) = 0;
v (end,:) = 1;
v (:, 1) = (0: h: 1) ';
v (:, end) = (0: h: 1) ';
pause (0.5);
surf (v);
i = i + 1;
end;
Результат роботи програми:
i = 1
d = 0.2250
d = 0.0875
d = 0.0500
d = 0.0344
d = 0.0297
d = 0.0245
d = 0.0199
d = 0.0175
d = 0.0154
d = 0.0137
d = 0.0120
d = 0.0108
d = 0.0093
HELM обчислювався методом МОНТЕ-КАРЛО (АЛГОРИТМ "Блукання по СІТЦІ") РІШЕННЯ задачі Діріхле для рівняння Гельмгольца в задану точку (x, y) Прямокутна область D, визначення меж.
Код програм:
ЗАПИСИ КРАЙОВИХ УМОВ У завдання
Діріхле для рівняння Гельмгольца
function yp = funch (x, y);
if x = 0, yp = y; end;
if y = 0, yp = 0; end;
if y = 1, yp = 1; end;
if x = 1, yp = y; end;
function [z1, z2, z3] = helm (c, fun, xm, ym, gr, x0, y0, h, n);
HELM обчислювався методом МОНТЕ-КАРЛО (АЛГОРИТМ "Блукання по СІТЦІ")
РІШЕННЯ задачі Діріхле для рівняння Гельмгольца в ЗАДАНОЇ
ТОЧЦІ (x, y) Прямокутна область D, визначення меж
0 <= x <= xm, 0 <= y <= ym
(УЧП) Uxx + Uyy-c * U = F (x, y)
(ГУ) U | г = G (x, y)
Вхідні дані:
c - КОЕФІЦІЄНТ (функціональний) ЛІВОЇ ЧАСТИНИ УЧП;
fun - ФУНКЦІЯ F (x, y) у правій частині УЧП (ФУНКЦІЯ КОРИСТУВАЧА);
xm, ym - КОРДОНУ Прямокутна область;
gr - Граничні умови (ФУНКЦІЯ КОРИСТУВАЧА);
x0, y0 - КООРДИНАТИ ТОЧКИ, В ЯКІЙ шукає рішення;
h - КРОК СІТКИ (задається користувачем);
n - ЧИСЛО Траєкторія.
Вихідні дані:
z1 - Наближені значення РІШЕННЯ;
z2 - ймовірність помилки;
z3 - ВЕРХНЯ КОРДОН ПОМИЛКИ.
Звернення: [z1, z2, z3] = helm (c, fun, xm, ym, gr, x0, y0, h, n)
global z
z = [];
i0 = round (x0 / h);
j0 = round (y0 / h);
im = round (xm / h);
jm = round (ym / h);
disp ('')
disp ('Зачекайте, йде розрахунок.')
for count = 1: n,
x = x0; y = y0;
i = i0; j = j0;
q = 1;
tmp = 4 + eval (c) * h ^ 2;
s = h ^ 2 * eval (fun) / tmp;
while all ([i, j, im-i, jm-j]),
p = [0,1 / 4]; p = [p, p (2)];
p = [p, 1 / 4]; p = [p, p (4)];
alf = rand;
pp = max (find (alf> cumsum (p)));
if pp == 1, j = j +1; end
if pp == 2, j = j-1; end
if pp == 3, i = i +1; end
if pp == 4, i = i-1; end
x = i * h; y = j * h;
q = q * 4/tmp;
s = s + q * h ^ 2 * eval (fun) / tmp;
end
s = s + q * feval (gr, x, y);
z = [z, s];
end
disp ('');
disp ('РІШЕННЯ ЗАВДАННЯ:');
disp ('=============================');
disp ('')
disp ('при числі траєкторій'); disp (n);
disp ('значення в точці з координатами');
disp ('x0 y0');
disp ([x0, y0]);
z1 = mean (z); disp ('Наближене рішення -'); disp (z1);
z2 = 0.6745 * std (z) / sqrt (n); disp ('ймовірні помилки -'); disp (z2);
z3 = z2 * 3/0.6745; disp ('ВЕРХНЬОЇ КОРДОНУ ПОМИЛКИ -'); disp (z3);
ЗВЕРНЕННЯ ДО ФУНКЦІЇ HELM:
global z
c = '1 ';
f = 0 ";
xm = 1; ym = 1;
gr = 'funch';
x0 = 0.6; y0 = 0.7;
h = 0.1;
n = 600;
[Z1, z2, z3] = helm (c, f, xm, ym, gr, x0, y0, h, n);
Результат роботи програми:
РІШЕННЯ ЗАДАЧІ:
при числі траєкторій 600
значення в точці з координатами x0 y0 0.6000 0.7000
Наближене рішення - 0.2958
Ймовірні помилки - 0.0089
ВЕРХНЬОЇ КОРДОНУ ПОМИЛКИ - 0.0397