Банк рефератов содержит более 364 тысяч рефератов, курсовых и дипломных работ, шпаргалок и докладов по различным дисциплинам: истории, психологии, экономике, менеджменту, философии, праву, экологии. А также изложения, сочинения по литературе, отчеты по практике, топики по английскому.
Полнотекстовый поиск
Всего работ:
364150
Теги названий
Разделы
Авиация и космонавтика (304)
Административное право (123)
Арбитражный процесс (23)
Архитектура (113)
Астрология (4)
Астрономия (4814)
Банковское дело (5227)
Безопасность жизнедеятельности (2616)
Биографии (3423)
Биология (4214)
Биология и химия (1518)
Биржевое дело (68)
Ботаника и сельское хоз-во (2836)
Бухгалтерский учет и аудит (8269)
Валютные отношения (50)
Ветеринария (50)
Военная кафедра (762)
ГДЗ (2)
География (5275)
Геодезия (30)
Геология (1222)
Геополитика (43)
Государство и право (20403)
Гражданское право и процесс (465)
Делопроизводство (19)
Деньги и кредит (108)
ЕГЭ (173)
Естествознание (96)
Журналистика (899)
ЗНО (54)
Зоология (34)
Издательское дело и полиграфия (476)
Инвестиции (106)
Иностранный язык (62792)
Информатика (3562)
Информатика, программирование (6444)
Исторические личности (2165)
История (21320)
История техники (766)
Кибернетика (64)
Коммуникации и связь (3145)
Компьютерные науки (60)
Косметология (17)
Краеведение и этнография (588)
Краткое содержание произведений (1000)
Криминалистика (106)
Криминология (48)
Криптология (3)
Кулинария (1167)
Культура и искусство (8485)
Культурология (537)
Литература : зарубежная (2044)
Литература и русский язык (11657)
Логика (532)
Логистика (21)
Маркетинг (7985)
Математика (3721)
Медицина, здоровье (10549)
Медицинские науки (88)
Международное публичное право (58)
Международное частное право (36)
Международные отношения (2257)
Менеджмент (12491)
Металлургия (91)
Москвоведение (797)
Музыка (1338)
Муниципальное право (24)
Налоги, налогообложение (214)
Наука и техника (1141)
Начертательная геометрия (3)
Оккультизм и уфология (8)
Остальные рефераты (21697)
Педагогика (7850)
Политология (3801)
Право (682)
Право, юриспруденция (2881)
Предпринимательство (475)
Прикладные науки (1)
Промышленность, производство (7100)
Психология (8694)
психология, педагогика (4121)
Радиоэлектроника (443)
Реклама (952)
Религия и мифология (2967)
Риторика (23)
Сексология (748)
Социология (4876)
Статистика (95)
Страхование (107)
Строительные науки (7)
Строительство (2004)
Схемотехника (15)
Таможенная система (663)
Теория государства и права (240)
Теория организации (39)
Теплотехника (25)
Технология (624)
Товароведение (16)
Транспорт (2652)
Трудовое право (136)
Туризм (90)
Уголовное право и процесс (406)
Управление (95)
Управленческие науки (24)
Физика (3463)
Физкультура и спорт (4482)
Философия (7216)
Финансовые науки (4592)
Финансы (5386)
Фотография (3)
Химия (2244)
Хозяйственное право (23)
Цифровые устройства (29)
Экологическое право (35)
Экология (4517)
Экономика (20645)
Экономико-математическое моделирование (666)
Экономическая география (119)
Экономическая теория (2573)
Этика (889)
Юриспруденция (288)
Языковедение (148)
Языкознание, филология (1140)

Лабораторная работа: Применение численных методов для решения уравнений с частными производными

Название: Применение численных методов для решения уравнений с частными производными
Раздел: Рефераты по математике
Тип: лабораторная работа Добавлен 11:45:52 16 августа 2010 Похожие работы
Просмотров: 97 Комментариев: 2 Оценило: 0 человек Средний балл: 0 Оценка: неизвестно     Скачать

САНКТ-ПЕТЕРБУРГСКИЙ УНИВЕРСИТЕТ ПУТЕЙ СООБЩЕНИЯ

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

ОТЧЕТ

ПО ВЫПОЛНЕННОЙ КУРСОВОЙ РАБОТЕ

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

«Применение численных методов для решения Уравнений с частными производными»

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

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.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) 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.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'

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.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(' 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.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(' 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.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

Оценить/Добавить комментарий
Имя
Оценка
Комментарии:
Где скачать еще рефератов? Здесь: letsdoit777.blogspot.com
Евгений07:56:43 19 марта 2016
Кто еще хочет зарабатывать от 9000 рублей в день "Чистых Денег"? Узнайте как: business1777.blogspot.com ! Cпециально для студентов!
08:20:31 29 ноября 2015

Работы, похожие на Лабораторная работа: Применение численных методов для решения уравнений с частными производными

Назад
Меню
Главная
Рефераты
Благодарности
Опрос
Станете ли вы заказывать работу за деньги, если не найдете ее в Интернете?

Да, в любом случае.
Да, но только в случае крайней необходимости.
Возможно, в зависимости от цены.
Нет, напишу его сам.
Нет, забью.



Результаты(150530)
Комментарии (1836)
Copyright © 2005-2016 BestReferat.ru bestreferat@mail.ru       реклама на сайте

Рейтинг@Mail.ru