САНКТ-ПЕТЕРБУРГСКИЙ УНИВЕРСИТЕТ ПУТЕЙСООБЩЕНИЯ
Кафедра «Прикладная математика»
ОТЧЕТ
ПО ВЫПОЛНЕННОЙ КУРСОВОЙ РАБОТЕ
Предмет «Численные методы»
«Применение численных методов длярешения Уравнений с частными производными»
Санкт-Петербург 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.00001.5000 2.0000 2.5000];
Y=[0.0378 0.0653 0.37891.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
P1 = -1.0362*X+2.5896
P2 = -2.3490*X2+7.1853*X+-4.4574
P3 = 2.8692*X3-15.2604*X2 + 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.00001.5000 2.0000 2.5000];
Y=[0.0378 0.0653 0.37891.0353 0.5172 0.9765];
cs = spline(X,[0 Y0]);
xx =linspace(0,2.5);
plot(X,Y,'*m',xx,ppval(cs,xx),'-k');
h=0.5
esstestvenniispline
A=[4 2 0 00 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 24]
B=[6*(Y(2)-Y(1))/h0 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'
otsutstvieuzla
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))/h0 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) Y S(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.82980 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 00 -3.9096
B1 = -1.1444 2.0466 5.82000 0 -3.9096
B1 = -1.1444 2.0466 5.82000.8298 0 -3.9096
B1 = -1.1444 2.0466 5.82000.8298 -0.3528 -3.9096
S1 =
0.2496
0.1008
1.3940
0.1433
-1.1372
4.0529
/>
Лабораторнаяработа N2 «Построение алгебраических многочленов наилучшегосреднеквадратичного приближения»
X=[0.0000 0.50001.0000 1.5000 2.0000 2.5000];
Y=[0.0378 0.06530.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'
X Y X^2 X*Y X^2*Y X^3 X^4
0 0.0378 0 0 0 0 0
0.5000 0.0653 0.2500 0.03260.0163 0.1250 0.0625
1.0000 0.3789 1.0000 0.37890.3789 1.0000 1.0000
1.5000 1.0353 2.2500 1.55302.3294 3.3750 5.0625
2.0000 0.5172 4.0000 1.03442.0688 8.0000 16.0000
2.5000 0.9765 6.2500 2.44136.1031 15.6250 39.0625
7.5000 3.0110 13.7500 5.440210.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=[nTABL(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.00000.5000 1.0000 1.5000 2.0000 2.5000];
Y=[0.03780.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=[nTABL(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 00 0
0.5000 0.9368 0.25000.4684 0.2342 0.1250 0.0625
1.0000 0.6846 1.00000.6846 0.6846 1.0000 1.0000
1.5000 0.3551 2.25000.5327 0.7990 3.3750 5.0625
2.0000 0.5962 4.00001.1924 2.3848 8.0000 16.0000
2.5000 0.3766 6.25000.9416 2.3539 15.6250 39.0625
7.5000 3.9122 13.75003.8196 6.4565 28.1250 61.1875
S2 =
6.0000 7.5000 13.7500
7.5000 13.750028.1250
13.7500 28.125061.1875
T2 =
3.9122
3.8196
6.4565
coef2 =
1.0178
-0.4243
0.0718
coef2 =
0.0718 -0.42431.0178
delt2 =
0.1199
delt2 =
0.0719
/>
Численныеметоды решения задачи Коши для обыкновенных дифференциальных уравнений
Эйлера явная
function dy=func(x,y)
dy=2*x*y
clear
X=[0.000000.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(' X Y 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.01h=0.001
0 1.0000 1.00001.0000 1.0000
0.1000 1.0101 1.00001.0090 1.0099
0.2000 1.0408 1.02001.0387 1.0406
0.3000 1.0942 1.06081.0907 1.0938
0.4000 1.1735 1.12441.1683 1.1730
0.5000 1.2840 1.21441.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(' X Y 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.01h=0.001
0 1.0000 1.00001.0001 1.0000
0.1000 1.0101 1.01011.0101 1.0101
0.2000 1.0408 1.04101.0408 1.0408
0.3000 1.0942 1.09471.0942 1.0942
0.4000 1.1735 1.17451.1735 1.1735
0.5000 1.2840 1.28581.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',[00.5],1)
TABL =
0 1.0000 1.0000 1.0000 1.0000
0.1000 1.0101 1.0204 1.01111.0102
0.2000 1.0408 1.0629 1.04301.0410
0.3000 1.0942 1.1308 1.09771.0945
0.4000 1.1735 1.2291 1.17871.1740
0.5000 1.2840 1.3657 1.29161.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([01],[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
[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
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);
ifx=0,yp=y;end;
ify=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
(УЧП) 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)));
ifpp==1,j=j+1;end
ifpp==2,j=j-1;end
ifpp==3,i=i+1;end
ifpp==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