Compare commits
No commits in common. "eb98e8630ff57d502dbb92ca8685d1ec816e5ae1" and "6c20a4c07c44b9649c0d519c7bb8910230fd5304" have entirely different histories.
eb98e8630f
...
6c20a4c07c
File diff suppressed because one or more lines are too long
Binary file not shown.
Binary file not shown.
File diff suppressed because it is too large
Load Diff
@ -1,161 +0,0 @@
|
||||
\documentclass{article}
|
||||
\usepackage[utf8]{inputenc}
|
||||
\usepackage{biblatex}
|
||||
\addbibresource{library.bib}
|
||||
\usepackage{listings}
|
||||
\usepackage{amssymb}
|
||||
\usepackage{comment}
|
||||
\usepackage{graphicx,amsmath}
|
||||
\newcommand{\norm}[1]{\left\lVert#1\right\rVert}
|
||||
\usepackage{hyperref}
|
||||
\hypersetup{
|
||||
colorlinks=true,
|
||||
linkcolor=blue,
|
||||
filecolor=magenta,
|
||||
urlcolor=cyan,
|
||||
pdftitle={Overleaf Example},
|
||||
pdfpagemode=FullScreen,
|
||||
}
|
||||
|
||||
\title{Numerical Methods: Lecture 4. Conditioning. Floating point arithmetic and stability. Systems of linear equations.}
|
||||
\author{Konstantin Tikhonov}
|
||||
|
||||
\begin{document}
|
||||
|
||||
\maketitle
|
||||
|
||||
\section{Suggested Reading}
|
||||
|
||||
\begin{itemize}
|
||||
\item Lectures 12-19, 20-23 of \cite{trefethen1997numerical}
|
||||
\item Lectures 6-7 of \cite{tyrtyshnikov2012brief}
|
||||
\end{itemize}
|
||||
|
||||
\section{Exercises}
|
||||
|
||||
Deadline: 18 Nov
|
||||
|
||||
\begin{enumerate}
|
||||
|
||||
\item (3) Propose a numerically stable way to compute the function $f(x,a)=\sqrt{x+a}-\sqrt{x}$ for positive $x,\;a$.
|
||||
|
||||
\item (2) Consider numerical evaluation $\mathcal{C}=\tan(10^{100})$ with the help of arbitrary-precision arithmetic module \lstinline{mpmath}, which can be called as follows:
|
||||
\lstset{language=Python}
|
||||
\lstset{frame=lines}
|
||||
% \lstset{label={lst:code_direct}}
|
||||
\lstset{basicstyle=\ttfamily}
|
||||
\begin{lstlisting}
|
||||
from mpmath import *
|
||||
mp.dps = 64 # precision (in decimal places)
|
||||
mp.pretty = True
|
||||
+pi
|
||||
\end{lstlisting}
|
||||
What is the relative condition number of evaluating $\mathcal{C}$ w.r.t the input number $10^{100}$? How many digits do you need to keep at intermediate steps to evaluate $\mathcal{C}$ with 7-digit accuracy?
|
||||
\begin{comment}
|
||||
\item (3) Check, that the following function
|
||||
\lstset{language=Python}
|
||||
\lstset{frame=lines}
|
||||
\lstset{label={lst:code_direct}}
|
||||
\lstset{basicstyle=\ttfamily}
|
||||
\begin{lstlisting}
|
||||
import math
|
||||
def round_to_n(x, n):
|
||||
if x == 0:
|
||||
return x
|
||||
else:
|
||||
return round(x, -int(math.floor(math.log10(abs(x)))) + (n - 1))
|
||||
\end{lstlisting}
|
||||
rounds $x$ to $n$ significant digits.
|
||||
A sample program to compute $\sum_{k=1}^{3000}k^{-2}\approx 1.6446$ via consequent summation with rounding of intermediate results to 4 digits looks as follows:
|
||||
\lstset{language=Python}
|
||||
\lstset{frame=lines}
|
||||
\lstset{label={lst:code_direct}}
|
||||
\lstset{basicstyle=\ttfamily}
|
||||
\begin{lstlisting}
|
||||
res = 0
|
||||
for k in range(1,3001):
|
||||
res = round_to_n(res+1/k**2, 4)
|
||||
\end{lstlisting}
|
||||
Despite the absence of subtractions (and related precision loss), this code allows to get only two significant digits. Explain, why this happens and propose a more accurate way to compute this sum (maintaining the restriction of keeping only 4 digits of intermediate result).
|
||||
\end{comment}
|
||||
\item (4) Implement the function \lstinline{solve_quad(b, c)}, receiving coefficients $b$ and $c$ of a quadratic polynomial $x^2 + b x + c$, and returning a pair of equation roots. Your function should always return two roots, even for a degenerate case (for example, a call \lstinline{solve_quad(-2, 1)} should result into \lstinline{(1, 1)}). Additionally, your function is expected to return complex roots.
|
||||
|
||||
After checking ensuring that your algorithm sort of works, try it on the following 5 tests. Make sure that all of them pass.
|
||||
\lstset{language=Python}
|
||||
\lstset{frame=lines}
|
||||
\lstset{label={lst:code_direct}}
|
||||
\lstset{basicstyle=\ttfamily}
|
||||
\begin{lstlisting}
|
||||
tests = [{'b': 4.0, 'c': 3.0},
|
||||
{'b': 2.0, 'c': 1.0},
|
||||
{'b': 0.5, 'c': 4.0},
|
||||
{'b': 1e10, 'c': 3.0},
|
||||
{'b': -1e10, 'c': 4.0}]
|
||||
\end{lstlisting}
|
||||
|
||||
\item (5) Consider the polynomial $$
|
||||
w(x)=\Pi_{r=1}^{20}(x-r)=\sum_{i=0}^{20} a_i x^i
|
||||
$$ and investigate the condition number of roots of this polynomial w.r.t the coefficients $a_i$. Perform the following experiment, using \texttt{numpy} root-finding algorithm. Randomly perturb $w(x)$ by replacing the coefficients $a_i\to n_i a_i$, where $n_i$ is drawn from a normal distribution of mean $1$ and variance $\exp(-10)$. Show the results of $100$ such experiments in a single plot, along with the
|
||||
roots of the unperturbed polynomial $w(x)$. Using one of the experiments, estimate the relative and absolute condition number of the problem of finding the roots of $w(x)$ w.r.t. polynomial coefficients.
|
||||
|
||||
\item (10)
|
||||
Consider the least squares problem $Ax\approx b$ at
|
||||
$$
|
||||
A = \begin{bmatrix}
|
||||
1 & 1\\
|
||||
1 & 1.00001\\
|
||||
1 & 1.00001
|
||||
\end{bmatrix},\quad b = \begin{bmatrix}
|
||||
2 \\
|
||||
0.00001 \\
|
||||
4.00001
|
||||
\end{bmatrix}.
|
||||
$$
|
||||
|
||||
\begin{itemize}
|
||||
\item
|
||||
Formally, solution is given by
|
||||
\begin{equation}
|
||||
\label{ex}
|
||||
x = ( A^T A )^{-1} A^T b.
|
||||
\end{equation}
|
||||
Using this equation, compute the solution analytically.
|
||||
\item Implement Eq. (\ref{ex}) in \lstinline{numpy} in single and double precision; compare the results to the analytical one.
|
||||
\item Instead of Eq. (\ref{ex}), implement SVD-based solution to least squares. Which approach is numerically more stable?
|
||||
\item Use \lstinline{np.linalg.lstsq} to solve the same equation. Which method does this function use?
|
||||
\item
|
||||
What are the four condition numbers of this problem, mentioned in Theorem 18.1 of Ref. \cite{trefethen1997numerical}? Give examples of perturbations $\delta b$ and $\delta A$ that approximately attain those condition numbers?
|
||||
\end{itemize}
|
||||
\item (7)
|
||||
Let $$A = \begin{bmatrix}
|
||||
\epsilon & 1 & 0\\
|
||||
1 & 1 & 1\\
|
||||
0 & 1 & 1
|
||||
\end{bmatrix}$$
|
||||
\begin{itemize}
|
||||
\item Find analytically LU decomposition with and without pivoting for the matrix $A$.
|
||||
\item Explain, why can the LU decomposition fail to approximate factors $L$ and $U$ for $|\epsilon|\ll 1$ in finite-precision arithmetic?
|
||||
\end{itemize}
|
||||
|
||||
\item (6) Consider computing the function $f(n, \alpha)$ defined by $f(0,\alpha)=\ln(1+1/\alpha)$ and recurrent relation
|
||||
\begin{equation}
|
||||
f(n,\alpha)=\frac{1}{n}-\alpha f(n-1,\alpha).
|
||||
\end{equation}
|
||||
Compute $f(20, 0.1)$ and $f(20, 10)$ in standard (double) precision. Now, do the same exercise in arbitrary
|
||||
precision arithmetic:
|
||||
\lstset{language=Python}
|
||||
\lstset{frame=lines}
|
||||
\lstset{label={lst:code_direct}}
|
||||
\lstset{basicstyle=\ttfamily}
|
||||
\begin{lstlisting}
|
||||
from mpmath import mp, mpf
|
||||
mp.dps = 64 # precision (in decimal places)
|
||||
f = mp.zeros(1, n)
|
||||
f[0] = mp.log(1+1/mpf(alpha))
|
||||
for i in range(1, n):
|
||||
f[i] = 1/mpf(i) - mpf(alpha)*f[i-1]
|
||||
\end{lstlisting}
|
||||
\end{enumerate}
|
||||
Plot the relative difference between exact and approximate results, in units of machine epsilon \texttt{np.finfo(float).eps} for $\alpha=0.1$ and $\alpha=10$ as function of $n$. How would you evaluate $f(30, 10)$ without relying on the arbitrary precision arithmetic?
|
||||
\printbibliography
|
||||
\end{document}
|
||||
@ -1,162 +0,0 @@
|
||||
\documentclass[prb, notitlepage, aps, 11pt]{revtex4-2}%
|
||||
\usepackage[utf8]{inputenc}
|
||||
\usepackage[T2A]{fontenc}
|
||||
\usepackage[english, russian]{babel}
|
||||
\usepackage{amsmath}
|
||||
\usepackage{enumitem}
|
||||
\usepackage{amsmath}
|
||||
\usepackage{delimset}
|
||||
\usepackage[pdftitle = a]{hyperref}
|
||||
\usepackage{datetime}
|
||||
\usepackage{minted}
|
||||
\usemintedstyle{friendly}
|
||||
\usepackage[a]{esvect}
|
||||
\hypersetup{
|
||||
colorlinks=true,
|
||||
linkcolor=blue,
|
||||
filecolor=magenta,
|
||||
urlcolor=cyan,
|
||||
pdfpagemode=FullScreen,
|
||||
}
|
||||
\usepackage{microtype}
|
||||
|
||||
\newcommand{\framesep}{0.6em}
|
||||
\BeforeBeginEnvironment{minted}{\vspace{-1.6em}}
|
||||
\AfterEndEnvironment{minted}{\vspace{-0.5em}}
|
||||
|
||||
\begin{document}
|
||||
|
||||
\begin{center}
|
||||
версия от \today \quad \currenttime
|
||||
\end{center}
|
||||
|
||||
\title{\texorpdfstring{
|
||||
Численные методы, осень 2022\\
|
||||
Задание 4 [Число обусловленности. Числа с плавающей точкой и вычислительная устойчивость] \\
|
||||
Всего баллов: 37 \ Срок сдачи: 18 ноября
|
||||
}{}
|
||||
}
|
||||
|
||||
\maketitle
|
||||
|
||||
\section*{Рекомендованная литература}
|
||||
|
||||
\begin{itemize}
|
||||
\item Лекции 12--19, 20--23 из \cite{trefethen1997numerical}
|
||||
\item Лекции 6--7 из \cite{tyrtyshnikov2012brief}
|
||||
\end{itemize}
|
||||
|
||||
\section*{Упражнения}
|
||||
|
||||
\begin{enumerate}
|
||||
|
||||
\item (3) Предложите вычислительно устойчивый способ вычислить функцию
|
||||
$f(x,a)=\sqrt{x+a}-\sqrt{x}$
|
||||
при положительных $x$ и $a$.
|
||||
|
||||
\item (2) Вычислите $\mathcal{C}=\tan(10^{100})$ с помощью модуля \mintinline{python}{mpmath}, предназначенного для арифметики произвольной точности. Пример использования:
|
||||
%
|
||||
\begin{minted}[frame = lines, framesep = \framesep]{python}
|
||||
from mpmath import *
|
||||
mp.dps = 64 # точность (число десятичных цифр)
|
||||
mp.pretty = True
|
||||
+pi # pi — переменная из mpmath
|
||||
\end{minted}
|
||||
%
|
||||
Чему равно относительное число обусловленности при вычислении $\mathcal{C} = \mathcal{C}(10^{100})$? Сколько цифр нужно хранить в памяти при промежуточных вычислениях, чтобы получить $\mathcal{C}$ с~точностью в 7 значащих цифр?
|
||||
|
||||
\item (4) Реализуйте функцию \mintinline{python}{solve_quad(b, c)}, возвращающую корни приведённого квадратного уравнения $x^2 + b x + c = 0$. Корни могут повторяться или быть комплексными.
|
||||
|
||||
Когда вам покажется, что функция работает, запустите её на следующих пяти тестах. Добейтесь того, чтобы она правильно работала на каждом из них.
|
||||
%
|
||||
\begin{minted}[frame = lines, framesep = \framesep]{python}
|
||||
tests = [{'b': 4.0, 'c': 3.0},
|
||||
{'b': 2.0, 'c': 1.0},
|
||||
{'b': 0.5, 'c': 4.0},
|
||||
{'b': 1e10, 'c': 3.0},
|
||||
{'b': -1e10, 'c': 4.0}]
|
||||
\end{minted}
|
||||
|
||||
\item (5) Рассмотрите многочлен
|
||||
$$
|
||||
w(x) = \prod_{r=1}^{20} (x-r) = \sum_{i=0}^{20} a_i x^i
|
||||
$$
|
||||
и исследуйте число обусловленности его корней, выступающих в роли функций от коэффициентов $a_i$. Проведите эксперимент: случайным образом измените коэффициенты и найдите новые корни с помощью алгоритма из \texttt{numpy}.
|
||||
Коэффициенты изменяйте по правилу $a_i \to n_i a_i$, где $n_i$ подчиняются нормальному распределению с математическим ожиданием, равным 1, и дисперсией, равной $\exp(-10)$. Проведите 100 таких экспериментов и изобразите результаты на одном графике вместе с корнями исходного многочлена.
|
||||
|
||||
Оцените по одному из экспериментов абсолютное и относительное число обусловленности корней многочлена как функций его коэффициентов.
|
||||
|
||||
\item (10)
|
||||
Рассмотрим задачу наименьших квадратов --- $Ax\approx b$:
|
||||
$$
|
||||
A = \begin{bmatrix}
|
||||
1 & 1\\
|
||||
1 & 1.00001\\
|
||||
1 & 1.00001
|
||||
\end{bmatrix},\quad b = \begin{bmatrix}
|
||||
2 \\
|
||||
0.00001 \\
|
||||
4.00001
|
||||
\end{bmatrix}
|
||||
$$
|
||||
|
||||
\begin{itemize}
|
||||
\item Формально решение можно найти как
|
||||
%
|
||||
\begin{equation}
|
||||
\label{ex}
|
||||
x = ( A^T A )^{-1} A^T b.
|
||||
\end{equation}
|
||||
%
|
||||
Вычислите его по этой формуле аналитически.
|
||||
|
||||
\item Вычислите (\ref{ex}) с помощью
|
||||
\mintinline{python}{numpy}, используя числа одинарной и двойной точности; сравните результат c аналитическим.
|
||||
|
||||
\item Помимо формулы (\ref{ex}), реализуйте решение, основанное на сингулярном разложении. Какой способ вычислительно более стабильный?
|
||||
|
||||
\item Решите эту же задачу с помощью \mintinline{python}{np.linalg.lstsq}. Какой алгоритм использует эта функция?
|
||||
|
||||
\item Какие четыре числа обусловленности, относящиеся к этой задаче, упоминаются в теореме 18.1 из~\cite{trefethen1997numerical}? (Возможно, их требуется вычислить --- прим. пер.).
|
||||
Приведите примеры таких $\delta b$ и $\delta A$, при которых приблизительно достигаются оценки на $\norm{\delta x}$, даваемые числами обусловленности.
|
||||
\end{itemize}
|
||||
|
||||
\item (7)
|
||||
Пусть
|
||||
$$
|
||||
A = \begin{bmatrix}
|
||||
\epsilon & 1 & 0\\
|
||||
1 & 1 & 1\\
|
||||
0 & 1 & 1
|
||||
\end{bmatrix}
|
||||
$$
|
||||
\begin{itemize}
|
||||
\item Аналитически найдите LU-разложение матрицы $A$ с применением выбора главного элемента и без него.
|
||||
\item Объясните, почему при $|\epsilon|\ll 1$ мы можем неправильно оценить множители $L$ и $U$ в арифметике конечной точности.
|
||||
\end{itemize}
|
||||
|
||||
\item (6) Пусть функция $f(n, \alpha)$ определена следующим образом:
|
||||
%
|
||||
\begin{align*}
|
||||
f(n,\alpha) &= \frac{1}{n} - \alpha f(n-1,\alpha) \\
|
||||
f(0,\alpha) &= \ln(1+1/\alpha)
|
||||
\end{align*}
|
||||
%
|
||||
Вычислите $f(20, 0.1)$ и $f(20, 10)$ с помощью арифметики обычной (двойной) точности. Теперь сделайте то же самое в арифметике произвольной точности:
|
||||
%
|
||||
\begin{minted}[frame = lines, framesep = \framesep]{python}
|
||||
from mpmath import mp, mpf
|
||||
mp.dps = 64 # precision (in decimal places)
|
||||
f = mp.zeros(1, n)
|
||||
f[0] = mp.log(1 + 1/mpf(alpha))
|
||||
for i in range(1, n):
|
||||
f[i] = 1/mpf(i) - mpf(alpha) * f[i-1]
|
||||
\end{minted}
|
||||
%
|
||||
Постройте в единицах машинного эпсилон график относительной разности между точными и приближёнными результатами как функции от $n$. Сделайте это при $\alpha=0.1$ и при $\alpha=10$. Машинный эпсилон можно получить как \mintinline{python}{np.finfo(float).eps}. \\
|
||||
Как бы вы стали вычислять $f(30, 10)$ без арифметики произвольной точности?
|
||||
|
||||
\end{enumerate}
|
||||
|
||||
\bibliography{library.bib}
|
||||
\end{document}
|
||||
@ -1,118 +0,0 @@
|
||||
@book{trefethen1997numerical,
|
||||
title={Numerical linear algebra},
|
||||
author={Trefethen, Lloyd N and Bau III, David},
|
||||
volume={50},
|
||||
year={1997},
|
||||
publisher={Siam}
|
||||
}
|
||||
@book{robert2013monte,
|
||||
title={Monte Carlo statistical methods},
|
||||
author={Robert, Christian and Casella, George},
|
||||
year={2013},
|
||||
publisher={Springer Science \& Business Media}
|
||||
}
|
||||
@book{murphy2012machine,
|
||||
title={Machine learning: a probabilistic perspective},
|
||||
author={Murphy, Kevin P},
|
||||
year={2012},
|
||||
publisher={MIT press}
|
||||
}
|
||||
@book{boyd2004convex,
|
||||
title={Convex optimization},
|
||||
author={Boyd, Stephen and Boyd, Stephen P and Vandenberghe, Lieven},
|
||||
year={2004},
|
||||
publisher={Cambridge university press}
|
||||
}
|
||||
@book{trefethen2019approximation,
|
||||
title={Approximation Theory and Approximation Practice, Extended Edition},
|
||||
author={Trefethen, Lloyd N},
|
||||
year={2019},
|
||||
publisher={SIAM}
|
||||
}
|
||||
@book{devroye:1986,
|
||||
author = {Devroye, Luc},
|
||||
date = {1986)},
|
||||
description = {Non-Uniform Random Variate Generation},
|
||||
keywords = {algorithms generation random simulation},
|
||||
location = {New York},
|
||||
publisher = {Springer-Verlag},
|
||||
title = {Non-Uniform Random Variate Generation(originally published with},
|
||||
year = 1986
|
||||
}
|
||||
@book{williams2006gaussian,
|
||||
title={Gaussian processes for machine learning},
|
||||
author={Williams, Christopher K and Rasmussen, Carl Edward},
|
||||
volume={2},
|
||||
number={3},
|
||||
year={2006},
|
||||
publisher={MIT press Cambridge, MA}
|
||||
}
|
||||
@book{hairer1993solving,
|
||||
title={Solving ordinary differential equations. 1, Nonstiff problems},
|
||||
author={Hairer, Ernst and N{\o}rsett, Syvert P and Wanner, Gerhard},
|
||||
year={1993},
|
||||
publisher={Springer-Vlg}
|
||||
}
|
||||
@book{hairer1993solving2,
|
||||
title={Solving ordinary differential equations. 2, Stiff and differential-algebraic problems},
|
||||
author={Hairer, Ernst and N{\o}rsett, Syvert P and Wanner, Gerhard},
|
||||
year={1993},
|
||||
publisher={Springer-Vlg}
|
||||
}
|
||||
@book{tyrtyshnikov2012brief,
|
||||
title={A brief introduction to numerical analysis},
|
||||
author={Tyrtyshnikov, Eugene E},
|
||||
year={2012},
|
||||
publisher={Springer Science \& Business Media}
|
||||
}
|
||||
@book{amosov2003,
|
||||
title={Numerical Methods for Engineers},
|
||||
author={Amosov, AA and Dubinskiy YuA and Kopchenova, NV},
|
||||
year={2003},
|
||||
publisher={Izdatelstvo MEI}
|
||||
}
|
||||
@article{arbenz2012lecture,
|
||||
title={Lecture notes on solving large scale eigenvalue problems},
|
||||
author={Arbenz, Peter and Kressner, Daniel and Z{\"u}rich, DME},
|
||||
journal={D-MATH, EHT Zurich},
|
||||
volume={2},
|
||||
year={2012}
|
||||
}
|
||||
@article{trefethen1996finite,
|
||||
title={Finite difference and spectral methods for ordinary and partial differential equations},
|
||||
author={Trefethen, Lloyd Nicholas},
|
||||
year={1996},
|
||||
publisher={Cornell University-Department of Computer Science and Center for Applied~…}
|
||||
}
|
||||
@book{boyd2018introduction,
|
||||
title={Introduction to applied linear algebra: vectors, matrices, and least squares},
|
||||
author={Boyd, Stephen and Vandenberghe, Lieven},
|
||||
year={2018},
|
||||
publisher={Cambridge university press}
|
||||
}
|
||||
@article{halko2011finding,
|
||||
title={Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions},
|
||||
author={Halko, Nathan and Martinsson, Per-Gunnar and Tropp, Joel A},
|
||||
journal={SIAM review},
|
||||
volume={53},
|
||||
number={2},
|
||||
pages={217--288},
|
||||
year={2011},
|
||||
publisher={SIAM}
|
||||
}
|
||||
@book{demmel1997applied,
|
||||
title={Applied numerical linear algebra},
|
||||
author={Demmel, James W},
|
||||
year={1997},
|
||||
publisher={SIAM}
|
||||
}
|
||||
@article{dahlquist198533,
|
||||
title={33 years of numerical instability, Part I},
|
||||
author={Dahlquist, Germund},
|
||||
journal={BIT Numerical Mathematics},
|
||||
volume={25},
|
||||
number={1},
|
||||
pages={188--204},
|
||||
year={1985},
|
||||
publisher={Springer}
|
||||
}
|
||||
Loading…
Reference in New Issue
Block a user