Решить задачу методом бисекции: Две лестницы одна в 20 м длины, другая в 30 м, поставлены поперёк улицы и опираются своими концами на противоположные дома. Определить с точностью до 1 см ширину улицы, если точка пересечения лестниц находится на высоте 8 м над землей.

 

Решение: Рассмотрим треугольники ABC и AFG, они являются подобными (по двум сторонам и углу между ними). Следовательно, выполняются  следующие соотношения:

  или 

Найдем BC по теореме Пифагора:

AB²=X²+BC²;

AB=20;

;

;

          (1);

Аналогичным образом рассмотрим треугольники CEA и CFG, они тоже подобные. Следовательно:

  или 

;

;

  (2);

Из (1) выражаем X₁:

(3);

Так как X=X₁+X₂, то X₂=X-X₁ следовательно подставляем в это выражение формулу (3) и получаем:

                               (4).

Решаем полученное уравнение методом бисекции (методом половинного деления).

Так как длина самой короткой лестницы составляет 20 метров, то можно предположить, что ширина улицы меньше этого значения. Поэтому сузим диапазон поиска корней до интервала от 2 до 18.

Пусть на отрезке (a, b) находится корень уравнения (4), если функция f(x) непрерывная, то возможны два варианта:

1)     f(a) > 0 и f(b) < 0 – функция убывает;

2)     f(a) < 0  и f (b) > 0 – функция возрастает.

Вычисляем необходимое число итерация для нахождения корня уравнения с заданной точностью:

==11 итераций (шагов).

 

Текст программы с комментариями:

% a - начало диапазона,

% b - конец диапазона,

% c - точность.

function y=bisek(a, b, c)

% i - счетчик количества итераций.

for i = 1:1:11

% t – середина отрезка (a, b).

t=(a+b)/2;

    if (urav(a)<0 & urav(t)>0)

        b=t;

    end

    if (urav(t)<0 & urav(b)>0)

        a=t;

    end

    if (urav(a)>0 & urav(t)<0)

        b=t;

    end

    if (urav(t)>0 & urav(b)<0)

        a=t;

    end

    plot (a, 0, 'r+');

    i=i+1

    y=t

end

 

% вычисляемое уравнение

function ur=urav(x)

ur = ((x-((8*x)/((400-(x^2))^(0.5))))*((900-(x^2))^(0.5)))-(8*x);

 

Вызов функции осуществляется следующей командой:

>>bisek (2, 18, 0.01)

Результат работы программы (шаг итерации и значение корня):

i = 1            y = 10

i = 2            y = 14

i = 3            y = 16

i = 4            y = 17

i = 5            y = 16.5000

i = 6            y = 16.2500

i = 7            y = 16.1250

i = 8            y = 16.1875

i = 9            y = 16.2188

i = 10          y = 16.2031

i = 11          y = 16.2109

ans = 16.2109

 

График функции:

На графике плюсами обозначено последовательное приближение к корню уравнения.

В результате получаем ответ на вопрос задачи – ширина улицы составляет 16 метров 21 см.

 

Интегрирование. Вычислить интеграл по формулам трапеции и Симпсона с точность (1/2)*10^-3:  .

Вычисление интеграла по формуле трапеций.

         Суть метода в следующем. Разобьём отрезок [a, b] на n равных частей длины h = (b-a)/n точками x_i = a + ih (I = 0, 1, …, n), x₀ = a, x_n = b. Искомый интеграл равен сумме интегралов. На каждом элементарном отрезке [x_i-1; x_i] заменим подынтегральную функцию линейно функцией вида

y = A(x – x_i-1) + B(x – x_i),  xÎ[x_i-1; x_i] (i=1,2,…,n).

Она в граничных точках x_i-1, x_i принимает те же значения, что и функция f(x). Её график является ломаной линией, начальная, угловые и конечная точки которой принадлежат также графику подынтегральной функции f(x). С увеличение n число общих точек будет расти, а в результате ломаная станет всё больше приближаться к кривой y=f(x).

Площадь трапеции можно вычислить по формуле:

;

Погрешность данного метода составляет:

,

где .

n – число интервалов:

;

.

 

Текст программы для вычисления интеграла по формуле трапеции:

function Str = trapec()

a = 0;  % Нижняя граница интегрирования

b = 1;  % Верхняя граница интегрирования

R = (1/2)*10^(-3);  % Точность

 

% Определяем количество интервалов.

x=1;   % Значение икса.

% MM - вторая производная.

MM = (-1/(4*x)^(3/2))*cos(x) - (1/(x^(1/2))*sin(x)) - x^(1/2)*cos(x);

n = ceil(((((b-a)^3)*abs(MM))/(12*R))^0.5);

 

Str=0;  % В этой переменной накапливаем результат

for Xi=a:((b-a)/(n+1)):b

    Xi_1 = Xi + ((b-a)/(n+1));

    Str = Str + ((funk(Xi) + funk(Xi_1))/2)*(Xi_1 - Xi);

end

 

% вычисляемая функция

function ur = funk(x)

ur=(x^(1/2))*(cos(x));

 

Результаты работы программы (номер шага и значение интеграла на этом шаге):

i = 1            Str = 0.0243

i = 2            Str = 0.0464

i = 3            Str = 0.0725

i = 4            Str = 0.1016

i = 5            Str = 0.1333

i = 6            Str = 0.1670

i = 7            Str = 0.2023

i = 8            Str = 0.2387

i = 9            Str = 0.2759

i = 10          Str = 0.3136

i = 11          Str = 0.3513

i = 12          Str = 0.3886

i = 13          Str = 0.4253

i = 14          Str = 0.4610

i = 15          Str = 0.4954

i = 16          Str = 0.5281

i = 17          Str = 0.5588

ans = 0.5588

В результате получаем значение интеграла с заданной точностью – 0,5588.

 

 

Вычисление интеграла по формуле Симпсона.

         Суть метода. Разобьем отрезок [a, b] на чётное число n = 2m равных частей длины h = (b-a)/(2m). Границы элементарных отрезков сузить x₀ = a, x₁, x₂, …, x_2m = b, так что x₁ = x₀+h, x₂ = x₁ + h = x₀ + 2h, …, x_2m = x_2m-1 + h = x₀ + 2mh. Кривую y = f(x) на отрезке [x_2i-2; x_2i] длины 2р аппроксимируем параболой

y = A(x – x_2i-2)(x – x_2i-1) + B(x – x_2i-2)(x – x_2i) + C(x – x_2i-1)(x – x_2i),     (5)

таким образом, площадь под кривой приближённо заменим площадью под этой параболой.

Найдем коэффициенты A, B и C. Для этого в формуле (5) последовательно положим x = x_2i = x_2i-2 + 2h, x = x_2i-1 = x_2i-2 + h и x = x_2i-2. Получим соответственно

,      ,   .              (6)

Площадь под параболой есть

.

Поэтому имеем . Используя (6), найдём

Откуда найдем .

Погрешность по данному методу составляет

,

где =

n – число интервалов.

.

Текст программы для вычисления интеграла по формуле Симпсона:

function Str = sims()

a = 0;  % Нижняя граница интегрирования

b = 1;  % Верхняя граница интегрирования

R = (1/2)*10^(-3);  % Точность

 

% Определяем количество интервалов.

x=1;

% MMMM - четвертая производная

MMMM = ((-15*cos(x))/(16*(x^(7/2)))) - ((3*sin(x)) / (2*(x^(5/2)))) + ((3*cos(x)) / (2*(x^(3/2)))) + ((x*sin(x)) / (x^(1/2)))+(x^(1/2))*cos(x);

n = ceil(((((b-a)^5)*MMMM)/(180*R))^(1/4));

h = (b-a)/n;

y0 = funk(a);

Y = 0;

for m=1:2:n

    m=m;

    y1 = 4 * funk(h*m);

    y2 = 2 * funk(h*(m+1));

Y=Y+(y1+y2);

end

 

Str=(h/3)*(Y+y0);

% вычисляемая функция

function ur = funk(x)

ur=(x^(1/2))*(cos(x));

 

n=6;

h = (b-a)/n=(1-0)/6=0.1667;

 

Результаты работы программы:

Y = 2.7015

Y = 6.4670

Y = 10.0029

Str = (h/3)*Y= (0.1667/3)*10.0029 = 0.5558.

 

 

Метод наименьших квадратов

Для функции cos(px) на отрезке [-1,1] найти среди многочленов третьей степени тот многочлен, который дает наилучшее приближение по методу наименьших квадратов, если использовать значения функции в точках:

x₀ = -1; x₁ = -0.5; x₂ = 0; x₃ = 0.5; x₄ = 1.

 

x	-1	-0.5	0	0.5	1
y	-1	0	1	0	-1

 

;

 

n = 5 – число точек;

m = 3 – степень полинома.

Решая полученную систему методом Крамера, находим коэффициенты:

a₀ = 0.65699;                  a₁ = 0;         a₂ = -1.71429;       a₃ = 0

Получаем полином следующего вида:

;

Графики исходного выражения и полученного полинома.

 

Среднеквадратичное отклонение:

=0,2285;

Решить систему уравнений Ax=b методом Гаусса.

        

 

Применим простейший вариант метода Гаусса – схему единственного деления. Следуя Гауссу, сначала преобразуем квадратную матрицу коэффициентов при неизвестных к верхней треугольной матрице, из которой затем найдем решение системы.

Схема единственного деления состоит в следующем. Вначале исколючим неизвестное x₁ из всех уравнений системы, кроме первого. Разделив первое уравнение на коэффициент а₁₁ при x₁. Далее, вычитая полученное уравнение, соответственное на a₂₁, a₃₁, a₄₁ из (2) (3) и (4) уравнений, получим систему трех уравнений с тремя неизвестными.

На следующем шаге проделываем аналогичные преобразования. Исключив теперь x₂. Объединяя полученные уравнения, придем к эквивалентной системе с треугольной матрицей коэффициентов.

 

Текст программы на MatLAB для решения данной системы уравнений:

A=[   6 4.98 -1.97 3.96 -3.95;

9.07 -0.99 3.96 -0.99 12.94;

2.99 3.91 2.01 -1.96 0.93;

2.94 -8.93 0 1.92 10.96];

for n=0:2

if(A(n+1,n+1)~=0)

for i=0:2-n

s=A(4-i,1+n)/A(1+n,1+n);

for k=1+n:5

A(4-i,k)=A(4-i,k)-A(1+n,k)*s;

end

end

end

end

 

x4=A(4,5)/A(4,4)

x3=(A(3,5)-A(3,4)*x4)/A(3,3)

x2=(A(2,5)-A(2,3)*x3-A(2,4)*x4)/A(2,2)

x1=(A(1,5)-A(1,2)*x2-A(1,3)*x3-A(1,4)*x4)/A(1,1)

 

Результаты работы программы:

x4 = 0.2562;         x3 = 1.7754;         x2 = -0.9839;        x1 = 0.5721

 

Решение системы уравнений методом Зейделя:

Ax = b;

                   , ;

;

Критерий окончания итераций:

;

e  - заданная точность.

Текст программы:

A=[ 2.99 3.91 2.01 -1.96 0.93;

       2.94 -8.93 0 1.92 10.96;

       9.07 -0.99 3.96 -0.99 12.94;

       6 4.98 -1.97 3.96 -3.95;];

 

x=[0 0 0 0];

 

q=0;

E=1*10^(-4);

for i=1:4

    s=0;

    for j=1:4

        if(i~=j)

            s=s+sqrt((A(i,j)/A(i,i))^2);

        end

        if(q<s)

            q=s;

        end

    end

end

s

if (s<1)   

    p=1;

    y=x;

while(p==1)

    x(1)=(-x(2)*A(1,2)-x(3)*A(1,3)-x(4)*A(1,4)+A(1,5))/A(1,1);

    x(2)=(-x(1)*A(2,1)-x(3)*A(2,3)-x(4)*A(2,4)+A(2,5))/A(2,2);

    x(3)=(-x(1)*A(3,1)-x(2)*A(3,2)-x(4)*A(3,4)+A(3,5))/A(3,3);

    x(4)=(-x(1)*A(4,1)-x(2)*A(4,2)-x(3)*A(4,3)+A(4,5))/A(4,4);

    s=0;

    for j=1:4

        if(s<sqrt((x(j)-y(j))^2))

            s=sqrt((x(j)-y(j))^2)

        end

    end

    if(s<=(1-q)*E/q)

        p=0;

    end

    y=x

end

x

end

 

Результат работы программы: