Применение численных методов для решения уравнений с частными производными (86148)

Посмотреть архив целиком

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

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













ОТЧЕТ

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

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

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









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

По данным таблицы запишем и решим нормальную систему МНК-метода:


Случайные файлы

Файл
115778.rtf
147754.rtf
15212-1.rtf
114854.rtf
110220.doc




Чтобы не видеть здесь видео-рекламу достаточно стать зарегистрированным пользователем.
Чтобы не видеть никакую рекламу на сайте, нужно стать VIP-пользователем.
Это можно сделать совершенно бесплатно. Читайте подробности тут.