lectures/Komp_obr_SFedU/04_Pract.tex
Edward Emelianov 3648795f0c add 4, 6-1
2021-12-10 01:32:42 +03:00

583 lines
15 KiB
TeX
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

\documentclass[a4paper,12pt]{extarticle}
\usepackage{/home/eddy/ed, verbatim}
\title{Практикум \No4: системы уравнений, интегралы, производные}
\author{}\date{}\nocolon
\long\def\task#1{\par\noindent\leavevmode\refstepcounter{sect}\llap{\textbf{\thesect}\;}\indent\textit{#1}\par}
\def\t#1{{\upshape\ttfamily #1}}
\begin{document}
\maketitle
\section{Системы уравнений}
%
%
%
\task{Решить систему уравнений}
$$\left\{
\begin{aligned}
-x_1+x_2+2x_3&=10;\\
3x_1-x_2+x_3&=-20;\\
-x_1+3x_2+4x_3&=40.
\end{aligned}
\right.
$$
Представим ее в виде $\B{Ax=b}$.
Инициализируем постоянные:
\begin{verbatim}
A=[-1 1 2; 3 -1 1; -1 3 4];
b=[10; -20; 40];
\end{verbatim}
Нам необходимо проверить на вырожденность матрицу~$\B A$:
\begin{verbatim}
det(A)
ans = 10.000
\end{verbatim}
Теперь решить данную систему можно несколькими способами.
\begin{enumerate}
\item Через обратную матрицу.
$$\B A^{-1}\B{Ax}=\B A^{-1}\B b,\quad \Arr \quad
\B x=\B A^{-1}\B b.$$
В Octave это примет вид:
\begin{verbatim}
x = inv(A)*b
x =
1.00000
19.00000
-4.00000
\end{verbatim}
Проверим решение:
\begin{verbatim}
A*x
ans =
10.000
-20.000
40.000
\end{verbatim}
\item Метод Гаусса. Приведем к верхней треугольной форме расширенную
матрицу~$(\B A:\B b)$:
\begin{verbatim}
rref([A b])
ans =
1.00000 0.00000 0.00000 1.00000
0.00000 1.00000 0.00000 19.00000
0.00000 0.00000 1.00000 -4.00000
\end{verbatim}
Слева мы получили единичную матрицу, что значительно упрощает вычисления.
Однако, если бы матрица не имела нулей в правом верхнем углу, мы все равно
могли бы найти корни системы (обратный ход метода Гаусса).
\item Автоматический метод Гаусса. В данном случае необходимо лишь воспользоваться
уже известным вам оператором <<левого>> (или обратного) деления:
\begin{verbatim}
x = A\b
x =
1.0000
19.0000
-4.0000
\end{verbatim}
Далее мы будем использовать именно этот способ решения линейных систем уравнений.
\end{enumerate}
%
%
%
\task{Решить систему уравнений, заданную вырожденной матрицей}
$$\left\{
\begin{aligned}
x_1+3x_2+7x_3&=5;\\
-x_1+4x_2+4x_3&=2;\\
x_1+10x_2+18x_3&=12.
\end{aligned}
\right.
$$
\begin{verbatim}
A = [ 1 3 7; -1 4 4; 1 10 18];
b = [5; 2; 12];
det(A)
ans = 0
\end{verbatim}
Так как определитель матрицы коэффициентов равен нулю, невозможно найти обратную
матрицу. Однако, можно воспользоваться способом решения через псевдообратную
матрицу:
\begin{verbatim}
x = pinv(A)*b
x =
0.38498
-0.11033
0.70657
% check
A*x
ans =
5.0000
2.0000
12.0000
\end{verbatim}
Однако, изменим вектор $\B b$:
\begin{verbatim}
b = [3;6;0];
x = pinv(A)*b
x =
-1.08920
1.25117
-0.52347
% check
A*x
ans =
-1.0000
4.0000
2.0000
\end{verbatim}
В этом случае решение не будет точным (точнее, оно вообще не является решением
данной системы). Проверим, возможно ли найти общее решение данной системы уравнения,
приведя к верхней треугольной форме расширенную матрицу $(\B A:\B b)$:
\begin{verbatim}
rref([A b])
ans =
1.00000 0.00000 2.28571 0.00000
0.00000 1.00000 1.57143 0.00000
0.00000 0.00000 0.00000 1.00000
\end{verbatim}
Последняя строка содержит ненулевой элемент лишь в столбце свободных членов, что
однозначно свидетельствует об отсутствии решений данной системы уравнений.
\section{Степенные уравнения}
%
%
%
\task{Найти решение уравнения $2x^2-4x+5=0$.}
Для этого необходимо инициализировать полином набором коэффициентов и найти корни
командой~\verb'roots'.
\begin{verbatim}
p = [2 -4 5];
x = roots(p)
x =
1.0000 + 1.2247i
1.0000 - 1.2247i
\end{verbatim}
Итак, корни нашего уравнения: $x=1\pm1.2247i$. Точность вычислений Octave можно
задать явно командой~\verb'format'. Для отображения результата в виде рациональных
дробей можно указать следующее.
\begin{verbatim}
format rat
x =
1 + 4801/3920i
1 - 4801/3920i
\end{verbatim}
Вернуться к прежнему виду результатов можно командой~\verb'format short'.
\task{Найти корни полинома $p(x)=x^4+2x^3-3x^2+4x+5$ и получить его график на отрезке $[-4, 2]$.}
\begin{verbatim}
p = [1 2 -3 4 5];
x = roots(p)
x =
-3.18248 + 0.00000i
0.95560 + 1.11480i
0.95560 - 1.11480i
-0.72873 + 0.00000i
x=[-4:.05:2]; y=polyval(p,x);
plot(x,y)
\end{verbatim}
Нарисуем ось Х:
\begin{verbatim}
hold on
plot([-4 2], [0 0],'k')
\end{verbatim}
Команда \verb'hold on' позволяет <<дорисовать>> что-либо на уже имеющемся
графике. Буква~\verb"'k'" в параметре означает рисование черным цветом. Отключить
вывод на один и тот же график можно командой~\verb'hold off'.
%
%
%
\task{Найти решение уравнения $y=x^3+x^2-3x-3$.}
Зададим функцию:
\begin{verbatim}
f = inline("x^3+x^2-3*x-3");
\end{verbatim}
Функция \t{fsolve} позволяет решать нелинейные уравнения, и ее можно применить в т.ч. к решению
степенных уравнений. Необходимо задать начальное приближение для поиска. Задавая разные значения,
получим разные корни:
\begin{verbatim}
fsolve (f, 1)
ans = 1.7321
fsolve (f, 0)
ans = -1
fsolve (f, -2)
ans = -1.7321
\end{verbatim}
Можем проверить корни:
\begin{verbatim}
p=[1 1 -3 -3]
p =
1 1 -3 -3
roots(p)
ans =
1.7321
-1.7321
-1.0000
\end{verbatim}
А теперь попробуем решить этим же методом систему уравнений:
$$\begin{cases}
\e^{-\e^{-(x+y)}} = y(1+x^2),\\
x\cos y + y\sin x = 1/2.
\end{cases}
$$
Для начала конвертируем их к виду $F(x)=0$:
$$\begin{cases}
\e^{-\e^{-(x+y)}} - y(1+x^2) = 0,\\
x\cos y + y\sin x - 1/2 = 0.
\end{cases}
$$
Запишем функцию, позволяющую вычислить обе компоненты:
\verbatiminput{Materials4Pract/04/F.m}
Теперь попробуем найти решение, начиная с $(0,0)$:
\begin{verbatim}
fsolve(@F, [0 0])
ans =
0.35325 0.60608
\end{verbatim}
\section{Численное интегрирование, дифференциальные уравнения}
%
%
%
\task{Найти интеграл $\Int_0^3 x(\sin\frac{1}{x})\sqrt{|1-x|}\,dx$.}
Для вычисления подынтегральной функции в каждом узле интегрирования, нам необходимо задать функцию
\t{i1.m}:
\verbatiminput{Materials4Pract/04/i1.m}
Для интегрирования с оптимальным расчетом квадратур можно использовать функцию \t{quad}:
\begin{verbatim}
[q, ier, nfun, err] = quad (@i1, 0, 3)
ABNORMAL RETURN FROM DQAGP
q = 1.9819
ier = 1
nfun = 5061
err = 0.00000011522
\end{verbatim}
\t{q}~-- результат интегрирования, \t{ier}~-- код ошибки интегрирования (при нормальной процедуре
равен 0), \t{nfun}~-- количество узлов интегрирования, \t{err}~-- оценка ошибки интегрирования.
Здесь и во многих других функциях первым аргументом является либо строка с именем функции, либо
ссылка на нее (как в данном случае), либо inline-функция.
Еще примеры интегрирования. Квадратурная формула Гаусса--Конрода:
\begin{verbatim}
f = inline ("x.^3");
quadgk (f, 0, 1)
ans = 0.25000
\end{verbatim}
Квадратура Кленшоу--Куртиса (и бесконечный предел интегрирования):
\begin{verbatim}
f = @(x) x.^3 .* exp (-x);
quadcc (f, 0, Inf)
ans = 6.0000
\end{verbatim}
Квадратура Симпсона:
\begin{verbatim}
f = inline ("x.^3");
quadv(f, 0, 1)
ans = 0.25000
\end{verbatim}
Автоматический выбор квадратуры:
\begin{verbatim}
integral(f, 0, 1)
ans = 0.25000
\end{verbatim}
%
%
%
\task{Найти интеграл $\Int_0^5 (x^4+2x^2-1)dx$}
Можно посчитать интеграл и другим способом, если задан полином: определим коэффициенты полинома,
вычислим новый полином, являющийся интегралом нашего, а затем, вычитая первообразные, найдем
искомый интеграл:
\begin{verbatim}
c = [1 0 2 0 -1];
i = polyint(c);
I = polyval(i, 5) - polyval(i, 0)
I = 703.33
\end{verbatim}
Аналогичным образом мы можем вычислять производные:
\begin{verbatim}
d = polyder(c);
polyval(d, [1:5])
ans =
8 40 120 272 520
\end{verbatim}
%
%
%
\task{Вычислить интеграл $\Int_0^1 dx\Int_{-1}^1 \cos(\pi xy)\sqrt{x|y|}\,dy$}
Для двухмерного интегрирования воспользуемся функцией \t{dblquad}
\begin{verbatim}
I = dblquad(@(x, y) cos (pi*x.*y) .* sqrt (x.*abs(y)), 0, 1, -1, 1)
I = 0.30892
% OR
I = quad2d(@(x, y) cos (pi*x.*y) .* sqrt (x.*abs(y)), 0, 1, -1, 1)
I = 0.30892
% OR
[I err] = integral2(@(x, y) cos (pi*x.*y) .* sqrt (x.*abs(y)), 0, 1, -1, 1)
I = 0.30892
err = 0.00000030870
\end{verbatim}
Тройные интегралы~--- \t{triplequad} или \t{integral3}.
%
%
%
\task{Решить дифференциальное уравнение $\dot{x}=-\e^t x^2$ при $x(0)=2$.}
Запишем функцию, вычисляющую $\dot{x}$:
\verbatiminput{Materials4Pract/04/ode1.m}
Заданим аргумент $t\in[0,5]$ как вектор в 50 экземпляров
\begin{verbatim}
t = linspace(0,5,50);
x = lsode(@ode1, 2, t);
plot(t,x)
\end{verbatim}
\t{lsode} решает простейшее уравнение $\frac{dy}{dx}=f(x,y)$ при начальных условиях $y(0)$ по
заданному вектору~$x$.
%
%
%
\task{Решить методом Рунге--Кутты дифференциальное уравнение ван~дер~Поля}
$y''+\mu(1-y^2)y'+y=0$, $\mu>0$.
Для начала перепишем это уравнение с заменой $y_1=y$, $y_2=y_1'$:
$y_2'=\mu(1-y_1^2)y_2-y_1$. Для простоты примем~$\mu=1$.
Введем функцию, описывающую наше уравнение (ее необходимо ввести как новый
m-файл и сохранить под именем \t{vdp1.m}):
\verbatiminput{Materials4Pract/04/vdp1.m}
Теперь найдем решение уравнения и отобразим графики функции~$y$ и ее первой
производной:
\begin{verbatim}
[t, y] = ode45(@vdp1, [0 20], [2; 0]);
plot(t, y(:,1), '-', t, y(:,2), '--')
\end{verbatim}
Функция \verb'ode45' в качестве первого параметра требует имя функции, в которой
описано дифференциальное уравнение; второй параметр~--- интервал, в котором
изменяется аргумент искомой функции; третий аргумент~--- начальные условия для
функции и ее производной.
Возвращаемое значение~$y$ содержит два столбца: в первом находится искомая
функция, а во втором~--- ее первая производная.
Итак, для численного решения дифференциального уравнения в Octave необходимо
сначала представить это уравнение в виде линейной системы
$$\left\{\begin{aligned}
y_1'&=f_1(x,y_1,\ldots,y_n),\\
y_2'&=f_2(x,y_1,\ldots,y_n),\\
\cdots\\
y_n'&=f_n(x,y_1,\ldots,y_n).
\end{aligned}
\right.
$$
Затем функции $f_1$, \ldots, $f_n$ следует определить как строки специальной
функции, которая будет играть роль первого параметра функции, решающей данное
уравнение.
\section{Численное дифференцирование}
%
%
%
\task{Для ряда данных вычислить производную и построить график функции и производной}
\begin{verbatim}
x = [0:0.01:10];
y = y=x.^2.*sin(x)+sin(x/11)-tan(x*222)/cos(x);
\end{verbatim}
Простейший способ найти производную~--- воспользоваться методом разделенных разностей. Функция
\t{diff} вычисляет разность $y(x_{i+1})-y(x_i)$. Производную $y'(x)$ мы можем рассчитать в нулевом
приближении либо как $\frac{y(x_{i+1})-y(x_i)}{x_{i+1}-x_i}$, либо как
$\frac{y(x_i)-y(x_{i-1})}{x_i-x_{i-1}}$.
Попробуем оба способа. Учитывая то, что мы имеем равномерно распределенный ряд, вычисления
упрощаются.
\begin{verbatim}
dy1=[0 diff(y)]/0.01;
dy2=[diff(y) 0]/0.01;
plot(x,[y;dy1;dy2])
\end{verbatim}
Благодаря гладкости функции и большому шагу, мы практически не видим разницы. Однако, если мы в
10~раз уменьшим шаг, сдвиг уже будет иметь значение.
Кстати, мы можем и простейшим образом (трапециями) вычислить интеграл:
\begin{verbatim}
iy = [cumsum(y)];
plot(x,[y;dy1;iy])
legend("F", "dF", "iF")
\end{verbatim}
Если добавить ось X (\t{plot(x,[y;dy1;iy], x, zeros(size(x)))}), поведение интегральной кривой
отлично отразится на оригинальной функции.
%
%
%
\task{Найти производную зашумленного ряда данных.}
(не удалять предыдущие данные!)
О функции \t{polyder} мы уже упоминали. Она отлично подходит для тех наборов данных, которые можно
аппроксимировать полиномом. Давайте повторим предыдущие вычисления \t{y}, но добавим шум в 10дБ:
\begin{verbatim}
yn = awgn(y, 10, "measured");
plot(x,[y;yn], x, zeros(size(x)))
\end{verbatim}
Естественно, функции \t{diff} и \t{cumsum} в данном случае будут давать ужасный результат:
\begin{verbatim}
plot(x,[yn;[0 diff(yn)]/0.01])
\end{verbatim}
Попробуем аппроксимировать нашу кривую полиномом десятой степени и сравнить на графике (а потом
сравним с оригиналом):
\begin{verbatim}
p=polyfit(x,yn, 10);
plot(x,[yn; polyval(p,x)])
plot(x,[y; polyval(p,x)])
\end{verbatim}
Естественно, в самом начале (в районе нуля) шумы настолько велики, что аппроксимация получается,
мягко говоря, не очень. Но это все равно лучше, чем начальный зашумленный ряд.
Теперь вычисляем производную и сравним с предыдущей.
\begin{verbatim}
dp = polyder(p);
dyp=polyval(dp, x);
plot(x,[dy1;dyp])
\end{verbatim}
И еще:
\begin{verbatim}
plot(x,[y;yn;dy1;dyp])
\end{verbatim}
Можно попробовать разные степени полинома для аппроксимации этой функции, сравнив результаты.
Еще одним вариантом вычисления производной является функция \t{gradient}. Здесь можно
<<автоматически>> учесть шаг:
\begin{verbatim}
plot(x,[y;dy1;gradient(y,0.01)])
\end{verbatim}
А в случае неравномерно распределенных данных, мы можем задать вектор \t{x}.
\begin{verbatim}
x = [0 0.1 0.5 1 1.1 1.5 7 7.1 7.2 7.5 10 10.5 12 15 20 20.1 25 45 47 56 100];
y = x.^2.*sin(x)+sin(x/11)-tan(x*222)/cos(x);
dyy = gradient(y, x);
x1 = [0:0.1:100];
y1 = x1.^2.*sin(x1)+sin(x1/11)-tan(x1*222)/cos(x1);
plot(x, dyy, x1, [0 diff(y1)])
% и сравним с ходом оригинальной функции
plot(x, [y; dyy], x1, [y1; [0 diff(y1)]])
% who is who
legend("bad", "dbad", "ori", "dori")
\end{verbatim}
\task{Вычислите вторую производную предыдущей функции}
Для этого можно воспользоваться функцией \t{del2} (дискретный Лапласиан):
\begin{verbatim}
plot(x,[y;del2(y)*1e4])
\end{verbatim}
Не забываем, что т.к. мы вычисляем вторую производную, то интервал необходимо возвести в квадрат!
Естественно, N-ю производную мы можем вычислить и многократным вызовом функции \t{diff}, если
данные распределены равномерно.
\section{Задания для самостоятельного выполнения}
\begin{enumerate}
\item
Решите систему уравнений
$$\left\{\begin{aligned}
x_1+2x_2+3x_3&=1;\\
2x_1-x_2+4x_3&=2;\\
x_1-3x_2+x_3&=3.
\end{aligned}\right.
$$
%(1,$\pi$,0)
\item Решите систему уравнений
$$\left\{\begin{aligned}
x_1+x_2/2+x_3/3&=1;\\
x_1/2+x_2/3+x_3/4&=0;\\
x_1/3+x_2/4+x_3/5&=0.
\end{aligned}\right.
$$
Обратите внимание, что определитель матрицы коэффициентов {\tt det(A) = 4.6296e-04}.
Такие системы называются плохо обусловленными. Их решения сильно осциллируют
при малейших изменениях коэффициентов матрицы.
\item
Решите уравнение $x^7-2x^5+3x^3-4x=0$.
%( $0$,
%$\pm1.2848$,
%$0.9304 \pm 0.8313i$,
%$-0.9304 \pm 0.8313i$).
\item
Решите систему уравнений
$$
\left\{\begin{aligned}
\e^{x+y} &= \sin x;\\
\cos x &= \ln y - 1.
\end{aligned}\right.
$$
\item
Вычислите $\Int_0^1 \ln (x+1)\sin x\, dx$.
\item
Вычислите $\Int_{-1}^2 dx\Int_{-\pi}^0 dy\Int_0^1\frac{ln(xyz)}{\cos(xy)}dz$.
\item
Найдите решение уравнения ван~дер~Поля при $\mu=5$.
% (1.2350)
\item
Постройте график решения задачи Коши методом Рунге--Кутты на интервале
$[0,1]$ для уравнения $y'=x^3\sin y+1$ при $y(0)=0$.
\item
Найдите решение системы уравнений:
$$
\left\{\begin{aligned}
2(x-4)^2 + 7(y-8)^2 &= z^2;\\
5(x-1)^2 +1 + 2z^2 &= 4(y+3)^2;\\
x^2 + y^2 + z^2 &= 0.\\
\end{aligned}\right.
$$
\item Вычислите производную и интеграл для ряда данных $y=y(t)$.
\begin{verbatim}
t = [0 1 3 5 9 10 11 15 20 21 23 25 50 52 57 59 60 73 94 96 99 100];
y = [41.6 -0.4 7.6 -25.8 5.3 23.1 636.7 -46.7 -3.7 -29.1 96.6 3.3 -9.4 56.7
17.5 -17.1 17.4 4.3 -0.3 12.3 85.9 44.2];
\end{verbatim}
\item Вычислить производную и интеграл для функции $\left(\frac{\sin x}{x}\right)^2$ с отношением
сигнал-шум 10дБ на
промежутке $[-10,10]$.
\end{enumerate}
\end{document}