1605 Вычислительная математика
 

Новости/Публикации
Результаты контрольных работ ВМ.2010г

Примеры программ
LU - разложение с выбором главного элемента
QR-разложение. Процедура Грамма-Шмита.
QR-разложение. Отражения Хаусхолдера.
QR-разложение. Вращения Гивенса.
Приведение матрицы к бидиагональной форме.
Приведение комплексной бидиагональной матрицы к матрице с действительными неотрицательными числами.
Аппроксимация полиномом Лагранжа (Matlab)
Интерполяция кубическими сплайнами
Метод Эйлера
Метод Рунге-Кутта
Лабораторная раблта №1. Аппроксимация методом наименьших квадратов
Лабораторная работа №2. Интегрирующее звено.

Материалы по другим дисциплинам
Официальный сайт каф. ИТ-6
Полезные материалы (в том числе по дипл. проектированию)
Цифровая обработка сигналов

Список литературы

	Деммель Дж. Вычислительная линейная алгебра. - М.Мир, 2001. Книга известного американского математика-вычислителя представляет собой учебник повышенного уровня по вычислительным методам линейной алгебры, рядом особенностей выделяющийся среди изданий этого типа
	Турчак Л.И. Основы численных методов. - Физматлит, 2004. Содержит основные сведения о численных методах, необходимые для первоначального знакомства с предметом. Излагаются основы численных методов для систем линейных и нелинейных уравнений, а также дифференциальных и интегральных уравнений. Имеется много задач, примеров и алгоритмов для облегчения понимания логической структуры рассматриваемых методов и их использования в расчетах на компьютерах.
	Н. Бахвалов, Н. Жидков, Г. Кобельков Численные методы. - Бином, 2008. Классический учебник по численным методам, переработанный с учетом современных тенденций в вычислительных методах. В данном издании устранены неточности и опечатки, имевшиеся в предыдущих изданиях, упрощены некоторые доказательства.

LU - разложение с выбором главного элемента

function [P,L,U]=lu_decompose(A)
sza = size(A) ;
n = sza(1) ;
P = eye(n) ;
for iNum=1:n-1
    % перестановка
    jjNum = iNum ;
    for jNum=iNum:n
        if abs(A(jNum,iNum))>abs(A(jjNum,iNum))
            jjNum = jNum ;
        end
    end
    % перестановка строк
    if jjNum~=iNum
        tmp = A(iNum,:) ;
        A(iNum,:) = A(jjNum,:) ;
        A(jjNum,:) = tmp ;

        tmp = P(iNum,:) ;
        P(iNum,:) = P(jjNum,:) ;
        P(jjNum,:) = tmp ;
    end
    % вычисление столбца L
    for jNum=iNum+1:n
        A(jNum,iNum) = A(jNum,iNum)/A(iNum,iNum) ;
    end
    % обновление матрицы A
    for ii=iNum+1:n
        for jj=iNum+1:n
            A(ii,jj) = A(ii,jj) - A(ii,iNum)*A(iNum,jj) ;
        end
    end
end
L = eye(n,n) ;
for ii=1:n
    for jj=1:ii-1
        L(ii,jj) = A(ii,jj) ;
    end
end
U = zeros(n,n) ;
for ii=1:n
    for jj=ii:n
        U(ii,jj) = A(ii,jj) ;
    end
end

QR-разложение. Процедура Грамма-Шмита.

function [q,r] = qrd(a)
mn=size(a) ;
m = mn(1) ;
n = mn(2) ;
q = zeros(m,n) ;
r = zeros(n,n) ;
for i=1:n
    q(:,i)=a(:,i) ;
    for j=1:i-1
        %r(j,i)=q(:,j)'*a(:,i) ;
        r(j,i)=q(:,j)'*q(:,i) ;
        q(:,i) = q(:,i) - r(j,i)*q(:,j) ;
    end
    r(i,i) = sqrt(sum(q(:,i)'*q(:,i))) ;
    q(:,i) = q(:,i)/r(i,i) ;
end

QR-разложение. Отражения Хаусхолдера.

function [Q,R] = QRH(A)
m = size(A,1) ;
n = size(A,2) ;
EQ = 1 ;
for i=1:min(m-1,n)
    u = house(A(i:m,i)) ;
    Ps = eye(m-i+1)-2*u*u' ;
    P = eye(m) ;
    P(i:m,i:m) = Ps ;
    %A(i:m,i:n) = Ps*A(i:m,i:n) ;
    A = P*A ;
    EQ = EQ*P;
end
Q=EQ(:,1:n) ;
R=A(1:n,:);

function u=house(x)
u = x ;
u(1) = x(1)+sign(x(1))*sqrt(x'*x) ;
u = u/sqrt(u'*u) ;

QR-разложение. Вращения Гивенса.

function [Q,R] = QRG(A)
m = size(A,1) ;
n = size(A,2) ;
Q = 1 ;
for jj=1:n
    for ii=m-1:-1:jj
        RN = givens_rot(A(ii,jj),A(ii+1,jj),ii,m) ;
        A = RN*A ;
        Q = Q*RN' ;
    end
end
R = A(1:n,:) ;
Q = Q(:,1:n) ;

function R=givens_rot(xi,xj,ii,m)
xabs = sqrt(xi^2+xj^2) ;
c = xi/xabs ;
s = -xj/xabs ;
R2 = [c -s;s c] ;
R = eye(m) ;
R(ii:ii+1,ii:ii+1) = R2 ;

Приведение матрицы к бидиагональной форме.

function [Q,V,D] = ddiag(A)
% bidiagonal decomposition Q*A*V = D
n = size(A,1) ;
Q = eye(n) ;
V = eye(n) ;
for k=1:n-1
    % clear under diag
    u = house(A(k:n,k)) ;
    Ps = eye(n-k+1) - 2*u*u' ;
    P = eye(n) ;
    P(k:n,k:n) = Ps ;
    A = P*A ;
    Q = P*Q ;
    % clear above diag
    u = house(A(k,k+1:n)')' ;
    Ps = eye(n-k) - 2*u'*u ;
    P = eye(n) ;
    P(k+1:n,k+1:n) = Ps ;
    A = A*P ;
    V = V*P ;
end
D = A ;

Приведение комплексной бидиагональной матрицы к матрице с действительными неотрицательными числами.

function [U,R,V] = rd(A)
% decompose complex bidiagonal matrix
% A = U*R*V, R = bidiagonal matrix with
% positive real values
[M,N] = size(A) ;
R = A ;

% initial U,V values
U = 1 ;
V = 1 ;

proc_req = 1 ;
while proc_req
    proc_req = 0 ;

    % left (by row) correction
    for m=1:M
        Um = eye(M) ;
        for n = 1:N
            if abs(imag(R(m,n)))>1e-10
                Um(m,m) = R(m,n)/abs(R(m,n)) ; % exp(j*phi)
                R(m,:) = R(m,:)/Um(m,m) ;
                proc_req = 1 ;
                break ;
            end
        end
        U = U*Um ;
    end

    % right (by columns) correction
    for n=1:N
        Vn = eye(N) ;
        for m = 1:M
            if abs(imag(R(m,n)))>1e-10
                Vn(n,n) = R(m,n)/abs(R(m,n)) ; % exp(j*phi)
                R(:,n) = R(:,n)/Vn(n,n) ;
                proc_req = 1 ;
                break ;
            end
        end
        V = V*Vn ;
    end
end
R = real(R) ;

Аппроксимация полиномом Лагранжа

%function mlagrange
clc ;
clear all ;
X = (0:1.1:4*pi)' ;
Y = cos(X).^2 ;
N = length(Y) ;
%Y = Y + randn(N,1)*0.01 ;

Xp = (0.5:0.1:4*pi-.5) ;
Np = length(Xp) ;
Yp = zeros(Np,1) ;
for k=1:Np
    x = Xp(k) ;
    for n=1:N
        pl = plagrange(X,n,x)*Y(n) ;
        Yp(k) = Yp(k) + pl ;
    end
end
plot(0:0.1:4*pi,cos(0:0.1:4*pi).^2,'k-.'); hold on; plot(X,Y,'r^') ;
plot(Xp,Yp,'b-')  ;
legend('Original function      ',...
       'Points of approximation',...
       'Lagrange approximation '...
   ) ;
title('Lagrange polynomial approximation') ;
grid on ;

function Ln = plagrange(X,n,x)
numr = 1 ;
denm = 1 ;
for k=1:length(X)
    if k~=n
        numr = numr*(x-X(k)) ;
        denm = denm*(X(n)-X(k)) ;
    end
end
Ln = numr/denm ;

Интерполяция кубическими сплайнами

clc ;
clear all ;
x = [0   1  2  3 3.5] ;
y = [-1  5 -2  1 -8] ;

free_spline = 1 ;
k1 = tan(-0.25*pi) ;
k2 = tan(0.4*pi) ;
% h vector
h = (x(2:end)-x(1:end-1)) ;

% number of splines
n = length(y) - 1 ;

% compute c_i: C*c = f ==> c = inv(C)*f ;
C = zeros(n+1,n+1) ;
if free_spline
    C(1,1) = 2 ;
else
    C(1,1) = -2/3*h(1) ;
    C(1,2) = -1/3*h(1) ;
end
for k=2:n
    C(k,k-1) = h(k-1) ;
    C(k,k)   = 2*(h(k-1)+h(k)) ;
    C(k,k+1)   = h(k) ;
end
if free_spline
    C(n+1,n+1) = 2 ;
else
    C(n+1,n) = 1/3*h(n) ;
    C(n+1,n+1) = 2/3*h(n) ;
end
% b vector
f = zeros(n+1,1) ;
if free_spline==0
    f(1) = k1-(y(2)-y(1))/h(1) ;
end
for k=2:n
    f(k) = ((y(k+1)-y(k))/h(k)-(y(k)-y(k-1))/h(k-1))*3 ;
end
if free_spline==0
    f(n+1) = k2 -(y(n+1)-y(n))/h(n) ;
end
c = inv(C)*f ;

a = y(1:n) ;

b = zeros(n,1) ;
for k=1:n
    b(k) = (y(k+1)-y(k))/h(k)-h(k)/3*(c(k+1)+2*c(k)) ;
end

d = zeros(n,1) ;
for k=1:n
    d(k) = (c(k+1)-c(k))/(3*h(k)) ;
end

% draw spline
hold off ;
plot(x,y,'r^') ; grid on ; hold on ;
for ns=1:n
    hx = 0.1 ;
    cx = x(ns):hx:x(ns+1) ;
    cy = zeros(length(cx),1) ;
    for k=1:length(cx)
        xc = hx*(k-1) ;
        cy(k) = a(ns)+b(ns)*xc + c(ns)*xc^2 + d(ns)*xc^3 ;
    end
    plot(cx,cy,'b-') ;
end

Метод Эйлера

function [x,y,err] = euler(vardefs,func,x0,y0,h,N,prec_y)
eval(vardefs) ;
y = zeros(1,N) ;
x = x0:h:x0+h*(N-1) ;
y(1) = y0 ;
xv = x0 ;
yv = y0 ;
err = 0 ;
for n=2:N
    yv = yv + h*eval(func) ;
    y(n) = yv ;
    xv = xv + h ;
    err = err + abs(eval(prec_y)-yv)^2 ;
end
err = err/(N-1) ;

Метод Рунге-Кутта

function [x,y,err]=runge4(vardefs,func,x0,y0,h,N,prec_y)
eval(vardefs) ;
y = zeros(1,N) ;
x = x0:h:x0+h*(N-1) ;
y(1) = y0 ;
err = 0 ;
for n=2:N

    xv = x(n-1) ;
    yv = y(n-1) ;
    k1 = eval(func) ;

    xv = x(n-1)+h/2 ;
    yv = y(n-1)+h/2*k1 ;
    k2 = eval(func) ;

    xv = x(n-1)+h/2 ;
    yv = y(n-1)+h/2*k2 ;
    k3 = eval(func) ;

    xv = x(n-1)+h ;
    yv = y(n-1)+h*k3 ;
    k4 = eval(func) ;

    y(n) = y(n-1) + h/6*(k1+2*k2+2*k3+k4) ;

    xv = x(n) ;
    yp = eval(prec_y) ;
    err = err + abs(yp-y(n))^2 ;
end
err = err/(N-1) ;

Лабораторная работа №1. Интегрирующее звено.

clc ;
clear all ;

subplot(2,1,1) ;

defvars = 'T=5*1' ; % R*C
eval(defvars) ;
y0 = 2 ;
t0 = 0 ;
C1 = (y0-cos(t0)/(T^2+1)-T/(T^2+1)*sin(t0))/exp(-t0/T) ;
defvars = [defvars sprintf(';C1=%f',C1)] ;
t = 0:0.01:50 ;
y = cos(t)/(T^2+1)+T/(T^2+1)*sin(t)+C1*exp(-t/T) ;
plot(t,y,'k-') ;
grid on ;

[t,y1,err1] = euler(defvars,'(cos(xv)-yv)/T',0,y0,0.5,100,...
    'cos(xv)/(T^2+1)+T/(T^2+1)*sin(xv)+C1*exp(-xv/T)') ;
hold on ;
plot(t,y1,'r-') ;

[t,y2,err2] = runge4(defvars,'(cos(xv)-yv)/T',0,y0,0.5,100,...
    'cos(xv)/(T^2+1)+T/(T^2+1)*sin(xv)+C1*exp(-xv/T)') ;
hold on ;
plot(t,y2,'m-') ;

%title('y''=e^{-y}sin(x), y(0)=0')

subplot(2,1,2) ;
bar([err1 err2],'r')
grid on ;
set(gca,'XTickLabel',['  euler    '; ...
                      'Runge-Kutta']) ;
title('error values') ;