Застосування чисельних методів для вирішення рівнянь з приватними похідними

[ виправити ] текст може містити помилки, будь ласка перевіряйте перш ніж використовувати.

скачати

Санкт-Петербурзький університет шляхів сполучення

Кафедра «Прикладна математика»

ЗВІТ

По виконанню курсових робіт

Предмет «Чисельні методи»

«Застосування чисельних методів для вирішення Рівнянь з приватними похідними»

Санкт-Петербург 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

Додати в блог або на сайт

Цей текст може містити помилки.

Математика | Лабораторна робота
99.6кб. | скачати


Схожі роботи:
Дослідження чисельних методів вирішення нелінійних рівнянь
Застосування похідної та інтеграла на вирішення рівнянь і нерівностей
Чисельні методи для вирішення нелінійних рівнянь
Метод Гаусса для вирішення систем лінійних рівнянь
Вивчення теореми Безу для вирішення рівнянь n-го ступеня при n2
Вивчення теореми Безу для вирішення рівнянь n го ступеня при n 2
Інтегрування і пониження порядку деяких диференціальних рівнянь з вищими похідними
Розробка програмного забезпечення для вирішення рівнянь з однією змінною методом Ньютона дотичних
Особливості математичних методів застосовуваних для вирішення економічних задач
© Усі права захищені
написати до нас