Учебное пособие 1856
.pdf%- delta - допустимое отклонение для р0
%- epsilon - допустимое отклонение для эначений
%функции у
%- maxl - максимальное число итераций
%Выход- р - приближение Стеффенсена к нулю
%- Q - матрица, содержащая последовательность
%Стеффенсена
%Инициализацияматрицы R
R=zeros(maxl,3); for k=l:maxl
for j=2:3
% ВычислениезнаменателявформулеметодаНьютона-Рафсона nrdenom=feval(df ,R(k, j-1);
% ВычислениеприближенийНьютона-Рафсона if nrdenom==0
'деление на нуль в методе Ньютона-Рафсона' break
else
R(k, j)=R(k,j-l)-feval(f, R(k,j-l))/nrdenom; end
%Вычисление знаменателя в ускоренной процессе
%Эйткена
aadenom=R(k,3)-2*R(k,2)+R(k,l);
% Вычисление приближений ускоренного процесса Эйткена if aadenom==0
'деление на нуль в ускоренном процессе Эйткена' break
else R(k+l,l)=R(k,l)-(R(k,2)-R(k,l))^2/aadenom;
end end
% Конецпрограммыприпоявленииделениянануль if (nrdenom==0)|(aadenom==0)
break end end
ПрограммаP8 (методМюллера). Нахождениекорняуравненияf(x) = 0 стремя заданнымиразличныминачальнымиприближениямиРо, P1 и Р2.
function [p,y,err]=muller(f,p0,p1,p2,delta, epsilon,maxl)
%Вход- f - функция, вводимаякак строка'f '
%- р0, p1 и р2 - начальные приближения
%- delta - допустимое отклонение для р0, р1 и р2
%- epsilon - допустимое отклонение для значений функции у
%- maxl - максимальное число итераций
%Выход- р- приближениеМюллеракнулюфункцииf
%- у - значение функции у = f (р)
%- err - ошибка приближения к р
%ИнициализацияматрицРиY
Р=[р0 р1 р2]; Y=feval(f,P);
%Вычисление а и b for k=l:maxl
h0=P(l)-P(3); hl=P:(2)-P(3); e0=Y(l)-Y(3); e1=Y(2)-YX3); c=Y(3);
denom=hl*hO^2-hO*hl^2; a=(e0*hl-el*h0)/denom; b=(el*h0^2-e0*hl^2)/denom;
% Подавлениелюбыхкомплексныхкорней if b^2-4*a*c > 0
disc=sqrt(b^2-4*a*c); else
disc=0; end
% Нахождение наименьшегокорня if b < О
disc=-disc; end z=-2*c/(b+disc); p=P(3)+z;
% Сортировка входных Р для поиска двух ближайших к р if abs(p-P(2))<abs(p-P(1))
Q=[P(2) P(l) P(3)]; P=Q;
Y=feval(f,P);
end
if abs(p-P(3))<abs(p-P(2)) R=[P(1)P(3) P(2)]; P=R;
Y=feval(f,P); end
%Замена входного Р ближайшим от р рn
Р(3)-р;
Y(3) = feval (f, Р(3)); Y==Y(3);
%Вычисление критерия останова
err=abs(z);
relerr=err/(abs(p)+delta);
if (err<delta)|(relerr<delta)|(abs(y)<epsilon) break
end end
Программа P9 (обратная подстановка). Решение верхней треугольной системы линейных уравнений АХ = В методом обратной подстановки. Обращаться к методу только в том случае, когда диагональные элементы не равны нулю. Сначала вычисляем xN — bN/aNN. а затем используем правило
|
bk |
−∑Nj=k+1(akj xj) |
|
Xk= |
|
|
для N=N-1,N-2,…,1 |
|
|
||
|
|
ak |
|
function X=backsub(A,B) |
|||
% Вход |
|
- А - верхняя треугольная матрица размера n x n |
%- В - матрица размера n x 1
%Выход - X - решение системы линейных уравнений АХ= В
%НаходимразмерматрицыВиинициализируемX n=length(B);
X=zeros(n,l);
X(n)=B(n)/A(n,n); for k=n-l:-l:l
X(k)=(B(k)-A(k,k+l:n)*X(k+l:n))/A(k,k); end
Программа P10 (построение верхней треугольной матрицы для применения метода обратной подстановки). Чтобы найти решение системы АХ = В, сначала приводим расширенную матрицу [А|В] к верхней треугольной форме, а затемвыполняемобратнуюподстановку.
function X = uptrbk(A,B)
%Вход - А - невырожденная матрица размера N х N
%- В - матрица размера N х 1
%Выход- X - матрицаразмераN х1, содержащаярешениеАХ=В
%ИнициализацияX ивременноесохранениематрицы
[N N]=size(A);
X=zeros(N,l);
C=zeros(l,N+l);
% Видрасширеннойматрицы:Aug= [A|B] Aug=[A В];
for p=l:N-l
% Частный выбор главного элемента для столбца р
[Y,j]=mах(abs(Aug(р:N,р)));
% Меняем местами строки р и j C=Aug(p,:); Aug(p,:)=Aug(j+p-l,:); Aug(j+p-l,:)=C; if Aug(p,p)==0
'А вырождена. Нет единственного решения' break
end
% Процессисключениядлястолбцар for k=p+l:N
m=Aug(k,p)/Aug(p,p);
Aug(k,p:N+l)=Aug(k,p:N+l)-m*Aug(p,p:N+l); end
end
%Обратная подстановка в [UIY] с использованием программы
%P3.1
X=backsub(Aug(l:N,l:N),Aug(l:N,N+l));
Программа P11 (РА = LU: разложение с выбором главного элемента). Построение решения линейной системы АХ = В, где А — невырожденнаяматрица.
function X = lufact(A,B)
%Вход - А - матрица размера NхN
%- В - матрица размера N х 1
%Выход - X - матрица размера N х 1, содержащая решение
%АХ = В
%Инициализация X, Y, временное сохранение матрицы С и
%строк заданной матрицы перестановок R
[N,N]=size(A);
X=zeros(N,l);
Y=zeros(N,l);
C=zeros(l,N);
R=1:N; for p=l:N-l
%Находим главный элемент строки для столбца р
[maxl,j]=max(abs(A(p:N,p)));
%Меняем местами строки р и j
С=А(р,:); A(p,:)=A(j+p-1,:); A(j+p-1,:)=C;
d=R(p); R(p)=R(j+p-l);
R(j+p-1)= d; if A(p,p)==0
'А вырождена. Нет единственного решения' break
end
%Вычисление множителя и размещение под диагональю
%матрицы А
for k=p+l:N
mult=A(k,p)/A(p,p); A(k,p) = mult; A(k,p+l:N)=A(k,p+l:N)-mult*A(p,p+l:N);
end
end
% Решение для Y Y(1)=B(R(1));
for k=2:N
Y(k)=B(R(k))-A(k, 1:k-1)*Y(1:k-1);
end
% РeшениедляX X(N)=Y(N)/A(N,N); for k=N-l:-l:l
X(k)=(Y(k)-A(k,k+l:N)*X(k-1:N))/A(k,k); end
Программа P12 (итерация Якоби). Решение линейной системы
АХ = В, начиная с исходного значения X = Р0, и генерирования последовательности {Pk}, которая сходится к решению. Эффективным условием для применения метода является то, что А — строго диагонально доминирующая матрица.
function X=jacobi(A,B,P,delta, maxl)
%Вход - А - невырожденная матрица размера N х N
%- В - матрица размера N х 1
%- Р - матрица размера N х 1, начальное предположение
%- delta - допустимоеотклонениедляР
%- maxl - максимальное число итераций
%Выход - X - матрица размера N х 1: приближение Якоби к решению
%АХ = В
N = length(В); for k=l: maxl for j=l:N
X(j)=(B(j)-A(j,[1:j-1, j+1:N])*P([1:j-1, j+1,:N]))/A(j, j); end
err=abs(norm(X'-P)); relerr=err/(norm(X)+eps); P=X';
if(err<delta)|(relerr<delta) break
end end X=X ' ;
Программа P13 (итерация Гаусса-Зейделя). Решение линейной системы АХ = В, начиная с исходного значения X = P0 методом генерирования последовательности {Pk}, которая сходится к решению. Эффективным условием для применения метода является то, что А — строго диагонально доминирующая матрица.
function X=gseid(A,B,P,delta, maxl)
%Вход - А - невырожденная матрица размера N х N
%- В - матрица размера N х 1
%- Р - матрица размера N х 1, начальное предположение
%- delta - допустимое отклонение для Р
%- maxl - максимльное число итераций
%Выход-X матрица размера Nх1: приближение
%Гаусса-Зейделя к решению АХ = В
N = length(B); for k=l:maxl for j=l:N
if j==l
X(1)=(B(N)-A(1,2:N)*P(2:N))/A(1,1); elseif j==N
X(N=(B(N)-A(N,1:N-1)*(X(1:N-1)))/A(N,N); else
%Х содержит k-e приближение и Pk-1
X(j)=(B(j)-A(j,1,j-1)*X(1:j-1)-A(j,j+1:N)*P(J+1:N))/A(j,j); end
end err=abs(norm(X’-P);
relerr=err/(norm(X)+eps); Р=Х';
if(err<delta)|(relerr<delta) break
end
end X=X';
Программа P14 (нелинейная итерация Зейделя). Решение нелинейной системы неподвижной точки X = G(X) с заданным начальным приближением P0 и генерированием последовательности {Pk}, которая сходится к решению Р.
function [P,iter] = seidel(G,P,delta, maxl)
%Вход- G - нелинейнаясистема, записаннаявформеМ-файлаG.m
%- Р - начальное приближение к решению
%- delta - грань ошибки
%- maxl - число итераций
%Выход- Р- приближениеЗейделякрешению
%- iter - число потребовавшихся итераций
N=length(P); for k=l:maxl
Х=Р;
% X k-oe приближениекрешению
for j=l:N A=feval('G',X);
% ВыводчленовX помереихвычисления
X(j)=A(j) end
err=abs(norm(X-P)); relerr-err/(norm(X)+eps);
Р=Х; iter=k;
if(err<delta)|(relerr<delta) break
end end
Программа P15 (метод Ныотона-Рафсона). Решение нелинейной системы F(X) — 0 с одним заданным приближением P0 и генерированием последовательности {Pk}, которая сходится к решению Р.
function [Р,iter,err]=newdim(F,JF,P,delta,epsilon,maxl)
%Вход - F - система, записанная в форме М-файла F.m %- JF - матрица Якоби F, записанная в форме М-файла JF.M
%- Р - начальное приближение к решению
%- delta - допустимое отклонение для Р
%- epsilon - допустимое отклонение для F(P)
%- maxl - максимальное число итераций
%Выход - Р - приближенное решение
%- iter - число потребовавшихся итераций
%- err - ошибка вычисления Р
Y=feval(F,P); for k=l:maxl
J=feval(JF,P);
Q=P-(J\Y')';
Z=feval(F,Q); err=norm(Q-P);
relerr=err/(norm(Q)+eps);
P=Q;
Y=Z;
iter=k;
if (err<delta)|(relerr<delta)|(abs(Y)<epsilon) break
end end
Программа P16 (приближение Лагранжа). Программа вычисляет полином Лагранжа P(x)= ∑kN=0 yk LN ,k (x) основанный на
N + 1 точке (xk,yk) для k = 0,1,..., N. function [C,L]=lagran(X,Y)
%Вход - X - вектор абсцисс
%- Y - вектор ординат
%Выход- С- матрицакоэффициентовинтерполирующего
%полинома Лагранжа
%- L - матрица коэффициентов полинома Лагранжа w=length(X);
n=w-l; L=zeros(w,w);
%ФормированиекоэффициентовполиномаЛагранжа
for k=l:n+l
V=l;
for j=l:n+l if k~=j
V=conv(V,poly(X(j)))/(X(k)-X(j)); end
end
L(k,:)=V;
end
%Определениекоэффициентовинтерполирующегополинома
%Лагранжа
C=Y*L;
Программа P17 (интерполяционный полином Ньютона). Для построения и вычисления полинома Ньютона степени <N, который
проходит через точки (xk,yk) = (xkf(xk)), k = 0,l,...,N:
Р(х) = d0,0 + d1,1(x - х0) + d2,2(x - х0)(х - x1)+ + • • • + dN,N(x - х0)(х -x1)---(x- XN-1),
function [C,D]=newpoly(X,Y)
%Вход- X векторабсцисс
%- Y вектор ординат
%ВыходСвекторкоэффициентовинтерполирующего
%полинома Ньютона
%- D таблица разностных отношений
n=length(X);
D=zeros(n,n);
D(:,1)=Y';
% Использование формулы для формирования таблицы for j=2:n
for k=j:n
D(k, j)=(D(k, j-1)-D(k-l, j-l))/(X(k)-X(k-j+1); end
end
%Определениекоэффициентовинтерполирующего
%полиномаНьютона
C=D(n,n);
for k=(n-l):-l:l C=conv(C,poly(X(k))); m=length(C); C(m)=C(m)+D(k,k;
end
Программа P18 (приближение Чебышева). Построение и вычисление интерполирующего полинома Чебышева степени N на интервале [-1; 1], где полином
P(x)= ∑Nj=0 cj T j (x) |
|
||
построенпоузлам |
(2k +1)π |
|
|
xk = cos( |
) . |
||
2N +2 |
function [C,X,Y]=cheby(fun,n,a,b)
%Входfun - функция, заданнаяв видестроки, которуюприближают
%- N - степень интерполирующего полинома Чебышева
%- а - левая крайняя точка
%- b - правая крайняя точка
%Выход- С- массив коэффициентовполинома
%- X - массив абсцисс
%- Y - массив ординат
if nargin==2, a=-l; b=l; end d=pi/(2*n+2); C=zeros(l,n+l);
for k=l:n+l X(k)=cos((2*k-l)*d);
end X=(b-a)*X/2+(a+b)/2; x=X;
Y=eval(fun); for k =l:n+l
z=(2*k-l)*d; for j=l:n+l
C(j)=C(j)+Y(k)*cos((j-l)*z); end
end C=2*C/(n+l); C(l)=C(l)/2;
Программа |
P19 (построение |
линии |
методом |
наименьших |
||
квадратов). |
Построение |
линии |
методом наименьших |
квадратов |
||
у = Ах + |
В, которая |
соответствует |
N данным |
точкам |
(x1; y1) , . . . , ( x N , y N ) .
function [A,B]=lsline(X,Y)
%Вход- X - векторабсциссразмера1хп
%- Y - вектор ординат размера 1хп
%Выход- А- коэффициентприхвуравненииАх+ В
%- В - постоянный коэффициент в уравнении Ах + В xmean=mean(X);
ymean=mean(Y); sumx2=(X-xmean)*(X-xmean)'; sumxy=(Y-ymean)*(X-xmean)'; A=sumxy/sumx2;
B=ymean-A*xmean;
Программа P20 (построение полинома методом наименьших квадратов). Построение полинома степени М методом наименьших квадратов РМ(Х) = c1+ c2x + с3х2 +...+ смхм-1 + см+1хм, который
N
подогнан по N точкам { ( xk ,U k )}k =1 .
function С = lspoly(X,Y,M)
%Вход- X - вектор абсцисс размера 1хn
%- Y - вектор ординат размера 1хn
%- М - степень полинома, построенного
%методом наименьших квадратов
%Выход-С - коэффициенты полинома n=length(X);
B=zeros(l:M+l);
F=zeros(n,M+l);
%Заполнение столбцов матрицы F степенями X
for k=l:M+l
F(:,k)=X>.-(k-l);
end ' . % Решение линейной системы
A=F'*F;
B=F'*Y';
С=А\В;
C=flipud(C);
Программа P21 (смыкающий кубический сплайн). Построение и вычисление интерполяционного смыкающего кубического сплайна
N
S(x) по N+ 1 заданнойточке { ( xk , yk )}k =0 function S=csfit(X,Y,dxO,dxn)
%Вход- X - векторабсциссразмераIxn
%- Y - вектор ординат размера Ixn
%- dx0 = S' (х0) - граничноеусловиенапервуюпроизводную
%- dxn = S' (xn) - граничное условие на первую производную
%Выход- S: строкиS - этокоэффициентывпорядкеубывания
%для интерполяционного кубического сплайна
N=length(X)-l;
H=diff(X);
D=diff(Y)./H;
A=H(2:N-1);
B=2*(H(1:N-1)+H(2:N));
C=H(2:N);
U=6*diff(D);
%Смыкающий сплайн с ограничениями в крайних %точках
В(1)=В(1)-Н(1)/2; U(l)=U(l)-3*(D(l)-dxO); B(N-1)=B(N-1)-H(N)/2; U(N-l)=U(N-i)-3*(dxn-D(N)); for k=2:N-l
temp=A(k-l)/B(k-l);
B(k)=B(k)-temp*C(k-l);
U(k)=U(k)-temp*U(k-l);
end M(N)=U(N-1)/B(N-l);
for k=N-2:-l:l
M(k+l)=(U(k)-C(k)*M(k+2))/B(k);
end
M(l)=3*(D(l)-dxO)/H(l)-M(2)/2; M(N-H)=3*(dxn-D(N))/H(N)-M(N)/2;
for k=0:N-l S(k+l,l)=(M(k+2)-M(k+1))/(6*H(k+l));
S(k+l,2)=M(k+l)/:2; S(k+l,3)=D(k+i)-ft(k+l)*(2*M(k+l)+M(k+2))/6; S(k+l,4)=Y(k+l);
end
Программа P22 (тригонометрические полиномы). Построение тригонометрического полинома степени М вида
P(x)= a20 +∑Mj=1 (aj cos( jx) +bj sin( jx)) ,
построенного по N равномерно отстоящим значениям Xk = - π +2π kk/N, k = 1,2…N. В программе допускается, что 2М + 1 ≤
N.
function [A,B]=tpcoeff(X,Y,M)
%Входвекторравноотстоящихабсцисснаинтервале [-p1,pl]
%- вектор ординат
%- степень тригонометрического полинома
%Выходвектор, содержащийкоэффициентыприcos(jx)
%- вектор, содержащий коэффициенты при sin(jx)
if M>maxl M=maxl;
end
A=zeros(l,M+l);
B=zeros(l,M+l);
Yends=(Y(l)+Y(N+l))/2;
Y(l)=Yends;
Y(N+l)=Yends;
A(l)=sum(Y); for j=l:M
A(j+l)=cos(j*X)*Y';
B(j+l)=sin(j*X)*Y';
end A=2*A/N; B=2*B/N; A(l)=A(l)/2;
Программа P23 (дифференцирование, использующее пределы).
Чтобы получить численное приближение f'(х), генерируем последовательность
F’(x) ≈Dk = |
f (x +10−k h) − f (x +10−k h) |
, |
k = 0,1,..., n |
|||
|
||||||
|
2(10−k h) |
|
|
|
||
до тех пор, пока |
не |
будет |
выполнено неравенство |
|||
|Dn+1—Dn-1|≤|Dn—Dn-1| |
или |
|Dn — |
Dn1l |
не будет меньше |
допустимого отклонения, что является попыткой найти наилучшее приближение f ‘(x) ≈Dn
function [L,n]=difflim(f.x.toler)
%Bход - f - функция, вводимая в виде строки 'f'
%- х - точка дифференцирования
%- toler - допустимое отклонение для ошибки
%Выход - L=[H' D' Е']:
N=length(X)-l; maxl=fix((N-l)/2);
%Н - вектор длин шагов
%D - вектор приближений производной
%Е - вектор граней ошибок
%- n - координата "наилучшего приближения"
maxl=15;
h=1;
H(l)=h;
D(l)=(feval(f,x+h)-feval(f,x-h))/(2*h); Е(1)=0;
R(l)=0; for n=l:2
h=h/10;
H(n+l)=h; D(n+l)=(feval(f,x+h)-feval(f,x-h))/(2*h); E(n+l)-abs(D(n+l)-D(n)); R(n+l)=2*E(n+l)*(abs(D(n+l))+abs(D(n))+eps);
end n=2;
while((E(n)>E(n+l))&(R(n)>toler))& n<maxl h-h/10;
H(n+2)=h; D(n+2)=(feval(f,x+h)-feval(f,x-h))/(2*h); E(n+2)=abs(D(n+2)-D(n+l)); R(n+2)=2*E(n+2)*(abs(D(n+2))+abs(D(n+l))+eps); n=n+l;
end n=length(D)-1; L=[H’ D’ E’];
Программа P24 (дифференцирование с использованием экстраполяции). Чтобы найти численное приближение f'(x), генерируем таблицу приближений D(j, k) для k ≤j и используем f‘(x) ≈D(n, n) в качестве окончательного ответа. Приближения
D(j, k) хранятся в нижней треугольной матрице, первый столбец которой равен
D(j,0)= |
f (x + 2− j h) − f (x −2− j h) |
|||
2− j+1 h |
|
|
||
|
|
|
|
|
и элементы строки j имеют вид |
|
|
||
D(j,k)=D(j,k-1)+ |
D( j,k −1) −D( j −1,k −1) |
для 1 ≤k ≤j. |
||
|
|
4k −1 |
|
|
function [D,err,relerr,n]=diffext(f,x,delta,toler)
%Вход - f - функция, вводимая как строка 'f '
%- delta - допустимое отклонение для ошибки
%- toler - допустимое отклонение для относительной ошибки
%Выход- D - матрицаприближенийпроизводных
%- err - грань ошибки
%- relerr - грань относительной ошибки
%- п - координата "наилучшего приближения"
егг=1; relerr=l; h=l;
j=1;
D(l,l)=(feval(f,x+h)-feval(f,x-h))/(2*h); while relerr>toler & err>delta &j<12
h=h/2;
D(j+l,l)=(feval(f,x+h)-feval(f,x-h))/(2*h); for k=l:j
D(j+l,k+l)=D(j+l,k)+(D(j+l,k)-D(j,k))/((4^k)-1);
end err=abs(D(j-1,j+l)-D(j,j));
relerr=2*err/(abs(D(j+l,j+l))+abs(D(j,j))+eps);
j=j+1;
end [n,n]=size(D);
Программа P25 (дифференцирование, использующее N + 1 узел).
Программа предназначена для численного приближения f'(x) полиномомНьютонаN-й степени
Р(х) = a0 + a,(x - x0) + a2(X--X0)(X - x1) +
+аз(х — x0)(X - x1)(x — x2) +...+aN(X — X 0 ) . . . ( X — XN-I)
ииспользует значение f ‘(xо) ≈ Р'(хо) в качестве ответа. Метод следует использовать для точки XQ. Для вычисления f'(xk) fa P'(xk) точки можно перегруппировать {Xk, X0, ..., Xk-1, Xk+1,..., XN}.
function [A,df]=diffnew(X,Y)
%Вход- X - векторабсциссразмера Ixn
%- У - вектор ординат размера Ixn
%Выход- А- векторразмераIxn, которыйсодержиткоэффициенты
%полинома Ньютона N-й степени
%- df - приближение производной
А=Y; N=length(X); for j=2:K
for k=N:-l:j A(k)=(A(k)-A(k-l))/(X(k)-X(k-j+l));
end end x0=X(l); df=A(2); prod=l;
nl=length(A)-l; for k=2:nl
prod=prod*(xO-X(k)); df=df+prod*A(k+l);
end
Программа P26 (составная формула трапеций). Предназначена для приближения интеграла
∫b |
f (x)dx ≈ |
h |
( f (a) − f (b)) + h∑kM=−11 f (xk ) |
||
|
|||||
a |
|
2 |
|
|
|
посредством подсчета |
f(x) |
в М + 1 |
равноотстоящей точке |
||
xk = a + kh, k = О, 1, 2,..., М. Заметим, что x0 = а и хM = b. |
|||||
function s=traprl(f,a,b,M) |
|
|
|
какстрока 'f' |
|
% Вход- f - подынтегральнаяфункция, вводимая |
%- а и b - верхний и нижний пределы интегрирования
%- М - число подынтервалов
%Выход- s - суммаформулытрапеций
h=(b-a)/M; s=0;
for k=l:(M-l) x=a+h*k; s=s+feval(f,x);
end s=b*(feval(f,a)+feval(f,b))/2+h*s;
Программа P27 (составная формула Симпсона). Предназначена дляприближенияинтеграла
∫b |
f (x)dx ≈ |
h |
( f (a) + f (b)) + |
2h |
∑kM=−11 |
f (x2k ) +∑kM=1 f (x2k −1) |
||
|
3 |
|
||||||
a |
3 |
|
|
|
|
|
||
посредством подсчета |
f(x) в 2М |
+ 1 |
равноотстоящей точке |
хk, = а+ kh, k = 0, 1, 2,..., 2М. Заметим, чтоX0= аиx2M= b
function s=simprl(f,a,b,M)
%Вход- f - подынтегральнаяфункция, вводимаякакстрока
%- а и b - верхний и нижний пределы интегрирования
%- М - число подынтервалов