Методы вычислительной математики

 

Задание 1

Решить заданную систему  линейных уравнений методом простых  итераций.

 ;

Перепишем в виде:

А=;   B=;  X=. 

 

 

Ax=b  (1) преобразуем равносильным образом к виду

X=Bx+c  (2)

1 и 2 имеют одно множество решений.

X⁽⁰⁾;

X⁽¹⁾=Bx⁽⁰⁾+c;

X⁽²⁾=Bx⁽¹⁾+c;

……………..

X⁽ⁿ⁺¹⁾=Bx⁽ⁿ⁾+c; — Итерационная последовательность (в скобках номер итерации);

 

=x^* — этот вектор (Х^*) и будет являться решением системы 2, а значит и системы 1.

=;

 

Достаточное условие сходимости, метода итерации.

Если норма ||B||<1, то итерационная последовательность будет сходиться к истинному решению со скоростью геометрической прогрессии, знаменатель которой q=||B||, причем условие будет выполняться при любом начальном приближении. Чаще всего метод простой итерации производят в форме Яковлева.

 

 

 

 

 

 

 

 

 

Используя программу находим решение системы линейных уравнений

 

 

 

Задание №2.

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

       Задача №01. 

Х⁴+5Х-3=0

         Решение. Для локализации корней  применим графический способ. Преобразуем  исходное уравнение к следующему  эквивалентному виду:

Х⁴ = 3 – 5Х

        Построив  графики функций, определяем, что  у решаемого уравнения имеется  два  корня, которые находятся в интервале -2 -1,5 и 0,5<x<1

           1) -2 -1,5

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

 Для функции f(x)= x⁴ + 5x -3 имеем:

Первая производная: f’(x)= 4x³ + 5;

Вторая производная: f’’(x)= 12x²;

f’(-2)=-32-5=-27;

f’’(-1,5)=-13.5-5=-18.5;

Производные не меняют знак внутри интервалов областей определения.

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

           f-2)f’’(-2)0

M₁=min|f’(x)|=18,5;

M₂=max|f’’(x)|=48.

Зададим точность вычисления 0.0001 Вычисления  будем проводить по формуле:

X0,…Xn+1=Xn- ;

          Вычисление с помощью программных  средств:

- Метод бисекций. Ответ: -1,87570190. Кол-во итераций-13.

- Метод Ньютона. Ответ: -1,87571935. Кол-во итераций-3.

2) 0,5<x<1

Найдем корни уравнения  методом простых итераций. Для  корректного использования данного  метода, необходимо записать функцию  в виде:

X=φ(x);

X=;

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

Max| φ’(x)|<1;

φ(x)= - ;

|φ(x)|=;

max|φ(x)|=0.8<1;

Зададим точность вычисления 0.0001.

Вычисление с помощью  программных средств:

- Метод простых итераций. Ответ: 0,57795841 Кол-во итераций-3

- Метод бисекций. Ответ: 0,57772827. Кол-во итераций-13.

- Метод Ньютона. Ответ: 0,57772408. Кол-во итераций-3.

 

Задача №02. 

2arctgx-=0

Решение. Для локализации  корней применим графический способ. Преобразуем исходное уравнение  к следующему эквивалентному виду:

2arctgx =

        Построив  графики функций, определяем, что  у решаемого уравнения имеется  2 корня, которые находится в интервале -1 -0,5 и 0,5<x<1.

1)  -1 -0,5

       Найдем  корни уравнения методом бисекции.

A=-1;   B=-0,5;

C==-0,75;

f(c)=f(-0,75)=-0,101817≠0;

f(a)=f(-1)=-1,070796;

f(a)*f(c)=f(-1)*f(-0,75)=0,109025<0; значит принимаем a=c=-0,75;

Действия продолжаются до тех пор пока f(с) не станет равным 0, или разница между a и b не станет меньше заданной погрешности.

Зададим точность вычисления 0.0001

          Вычисление с помощью программных  средств:

- Метод бисекций. Ответ: -0,73371887. Кол-во итераций -14.

2)  0,5 1

       Найдем  корни уравнения методом бисекции.

A=1;   B=0,5;

C==0,75;

f(c)=f(0,75)=0,101817≠0;

f(a)=f(1)=1.070796;

f(a)*f(c)=f(1)*f(0,75)=0,109025<0; значит принимаем a=c=0,75;

Действия продолжаются до тех пор пока f(с) не станет равным 0, или разница между a и b не станет меньше заданной погрешности.

Зададим точность вычисления 0.0001

          Вычисление с помощью программных  средств:

- Метод бисекций. Ответ: 0,73371887. Кол-во итераций -14.

 

 

Задание №3.

       Используя обобщенные формулы трапеций и Симпсона вычислить определенные интегралы с заданной точностью. Проверку достижения требуемой точности проводить по правилу Рунге..

       Задача №01. (Метод Трапеций)

 dx

f(x) - непрерывная функция на отрезке

x₀  = a -нижний предел интегрирования.

x_n = b – верхний предел интегрирования.

X₀=0;

X_n=1;

Предположим что Y=h(₊f₁+f₂+…+f_n-1) –верное решение интеграла.

       Разобьем  отрезок на 5 частей  h=(1-0)/5=0,2 (n=5-кол-во разбиений)

f(x)=

f₀=f(X₀)=f(0)=0;

f₁=f(X₁)=f(0,2)=0,039479;

f₂=f(X₂)=f(0,4)=0,152203;

f₃=f(X₃)=f(0,6)=0,324252;

f₄=f(X₄)=f(0,8)=0,539193;

f₅=f(X₅)=f(1)=0,785398;

Y=0.2(+0,039479+0,152203+0,324252+0,539193+0,785398)=0.289685

 

 

 

 

       Задача №01. (Метод Симпсона)

 dx

     Метод Симпсона- график подынтегральной функции  приближается не ломанной линией, а маленькими параболами.

Предположим Y=h/3(f₀+f_2m+4_2i-1+2_2i)  верное решение интеграла.

Разобьем отрезок на 8 частей с шагом h=(1-0)/8=0,125;

f0=f(0)=0;

f1=f(0,125)=0,015544;

f2=f(0,25)=0,061245;

f3=f(0,375)=0,134539;

f4=f(0,5)=0,231824;

f5=f(0,625)=0,349125;

f6=f(0,75)=0,482626;

f7=f(0,875)=0,628976;

f8=f(1)=0,785398;

Y=0,125/3(0+0,785398+4(0,015544+0,134539+0,349125+0,628976)+2(0,06124++0,231824+0,482626))=0.285397

Проверка Аналитическим  методом

= -)dx=

=* – +)dx= – +)= +(+)| ₌ –

 

 

 

 

Задание №3.

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

       Задача №02.

, xϵ[0;1], h=0,1; N=8.

Решение: Используем программные  средства. Составим таблицу значений функции.

 

 

Задание №4.

     Найти приближенное  решение задачи Коши методом  Эйлера и Рунге-Кутта 4 порядка  на заданном отрезке с шагом  h=0.1.     

       Задача №01.

y’=;   y(2)=1; xϵ[2;2,25];

Задача Коши состоит в  следующем:

Найти решение y=y(x) системы или векторного уравнения, удовлетворяющее заданному начальному условию y(x₀)=y₀.

Расчетные формулы:

Метод Эйлера

 

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

Y_i+1=Y_i+(k₁+2k₂+2k₃+k₄);

k₁=f(X_i;Y_i);

k₂=f(X_i+;Y_i+);

k₃=f(X_i+;Y_i+);

k₄=f(X_i+h; Y_i+hk₃).

Метод Эйлера-Коши более  прост и требует меньшего количества вычислений, чем метод Рунге-Кутта, однако он менее точен; что является основным его недостатком.

Используя программные средства решаем уравнение.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

       Задача №02.

Найти приближенное решение задачи Коши y’’=f(x;y;y’), y(x₀)=y₀, y’(x₀)=y₀’или y’’’=f(x;y;y’;y’’), y(x₀)=y₀; y’(x₀)=y₀’; y’’(x₀)=y₀’’ методами Эйлера и Рунге-Кутта 4 порядка на отрезке [0;1] с шагом h=0,1 (h=0,01).

y’’+2y’+5y=-8sin2x; y(0)=2, y’(0)=6.

Для удобности запишем  уравнение в виде:

 

Решение: Используем программные  средства. Составим таблицу значений функции.

При уменьшении шага, h, точность результата вычисления увеличивается

 

program RungeKutt4;

uses Crt;

const

     NMax     =1;

     Space    ='               ';

type

       Mas      =array[1..NMax] of real;

var

     Y, Yo: mas;

     X, Xo, Xn, h:               real;

     i, j,  M :               integer;

     ch:                         char;

     procedure p(x:real; y:mas; var f: mas);

begin

     f[1]:=(y[1]*y[1]-x*x)/(2*x*y[1]);

end;

 

procedure Rgk(x,h:real; n:integer; var y:mas);

var

     k1, k2, k3, k4, y1:          mas;

     i:                       integer;

begin

     p(x, y, k1);

     for i:=1 to n do

        y1[i]:=y[i]+h*k1[i]*0.5;

     p(x+h*0.5, y1, k2);

     for i:=1 to n do

        y1[i]:=y[i]+h*k2[i]*0.5;

     p(x+h*0.5, y1,k3);

     for i:=1 to n do

     y1[i]:=y[i]+h*k3[i];

     p(x+h, y1, k4);

     for i:=1 to n do

     y[i]:=y[i]+h*0.166667*(k1[i]+2*k2[i]+2*k3[i]+k4[i]);

end;

 

begin

     ClrScr;

     writeln(Space, 'Приближенное решение задачи Коши');

     writeln(Space, '--------------------------------');

     writeln(Space, '  Метод  Рунге-Кутта 4 порядка   ');

     writeln(Space, '  ---------------------------   ');

     writeln                                           ;

     write  (       'Введите начальное значение Xo:  '); readln (Xo)  ;

     write  (       'Введите конец интервала Xn:     '); readln (Xn)  ;

     write  (       'Введите шаг h:                  '); readln (h)   ;

     writeln(       'Введите начальные условия:      ');

     for i:=1 to NMax do

begin

     write  (       ' Yo[',i,']=                     '); readln(Yo[i]);

end;

     M:= Trunc ((Xn-Xo)/h);

     X:= Xo;

     for i:=1 to NMax do Y[i]:=Yo[i];

     writeln;

     writeln (      'Результаты :                    ');

     for i:=1 to M do

begin

     Rgk (X, h, NMax, Y);

     X := X+h;

     write   (      ' X =', X:10:7);

     for j:=1 to NMax do

     write   (      'Y[',j,'] =' , Y[j]:10:7);

     writeln;

   end;

     writeln (      'Нажмите любую клавишу ...       ');

     ch := readkey;

end.

      

program RungeKutt4;

uses Crt;

const

     NMax     =2;

     Space    ='               ';

type

     Mas      =array[1..NMax] of real;

var

     Y, Yo: mas;

     X, Xo, Xn, h:               real;

     i, j,  M :               integer;

     ch:                         char;

     procedure p(x:real; y:mas; var f: mas);

begin

     f[1]:=y[2];

     f[2]:=-2*y[2]-5*y[1]-8*exp(-x)*sin(2*x);

end;

 

procedure Rgk(x,h:real; n:integer; var y:mas);

var

     k1, k2, k3, k4, y1:          mas;

     i:                       integer;

begin

     p(x, y, k1);

     for i:=1 to n do

        y1[i]:=y[i]+h*k1[i]*0.5;

     p(x+h*0.5, y1, k2);

     for i:=1 to n do

        y1[i]:=y[i]+h*k2[i]*0.5;

     p(x+h*0.5, y1,k3);

     for i:=1 to n do

     y1[i]:=y[i]+h*k3[i];

     p(x+h, y1, k4);

     for i:=1 to n do

     y[i]:=y[i]+h*0.166667*(k1[i]+2*k2[i]+2*k3[i]+k4[i]);

end;

 

begin

     ClrScr;

     writeln(Space, 'Приближенное решение задачи Коши');

     writeln(Space, '--------------------------------');

     writeln(Space, '  Метод  Рунге-Кутта 4 порядка   ');

     writeln(Space, '  ---------------------------   ');

     writeln                                           ;

     write  (       'Введите начальное значение Xo:  '); readln (Xo)  ;

     write  (       'Введите конец интервала Xn:     '); readln (Xn)  ;

     write  (       'Введите шаг h:                  '); readln (h)   ;

     writeln(       'Введите начальные условия:      ');

     for i:=1 to NMax do

begin

     write  (       ' Yo[',i,']=                     '); readln(Yo[i]);

end;

     M:= Trunc ((Xn-Xo)/h);

     X:= Xo;

     for i:=1 to NMax do Y[i]:=Yo[i];

     writeln;

     writeln (      'Результаты :                    ');

     for i:=1 to M do

begin

     Rgk (X, h, NMax, Y);

     X := X+h;

     write   (      ' X =', X:10:7);

     for j:=1 to NMax do

     write   (      'Y[',j,'] =' , Y[j]:10:7);

     writeln;

   end;

     writeln (      'Нажмите любую клавишу ...       ');

     ch := readkey;

end.

         

 

{Метод Ньютона}

program new;

uses CRT;

  const

    eps0=0.0001; m1=8,5; m2=48; x0=-2;

  var

    k:integer; x,eps:real;

  function f(x:real):real;

    begin

      f:=x*x*x*x+5*x-3

    end;

  function f1(x:real):real;

    begin

       f1:=10-exp(x)

     end;

  procedure Newton (eps:real; var x:real; var k:integer);

     var y:real;

     begin

        k:=0;

        repeat

           y:=x;

           x:=y-f(y)/f1(y);

           k:=k+1

        until abs (y-x) <eps

       end;

      begin

      x:=x0;

      eps:=sqrt(2*m1*eps0/m2);

      Newton(eps,x,k);

      writeln('количество итераций=',k:3);

      writeln('корень уравнения=',x:12:8)

end.

        

 

{Метод простых итераций}

program itera;

uses crt;

   const

      eps=0.001; x0=0,5;

   var

      x:real; k:integer; {ch:char;}

   function f(x:real):real;

      begin

            f:=(3-x*x*x*x)/5

      end;

   procedure siter (eps:real; var k:integer; var x:real);

   var       y:real;

   begin

             k:=0;

             y:=f(x);

   while abs(y-x)>=eps do

       begin

         k:=k+1;

         x:=y;

         y:=f(x)

       end

   end;

begin

  x:=x0;

  siter(eps,k,x);

writeln('кол-во итераций=',k:3);

writeln('корень уравнения=',x:12:8)

  {ch:=readkey;}

end.

 

{Метод бисекций или  деления отрезка пополам}

program popolam;

uses crt;

    const

      eps=0.0001; a=-2; b=-1,5;

    var

      x:real; k:integer;

    function f(x:real):real;

begin

      f:=x*x*x*x+5*x-3

      end;

procedure bisect (a,b,eps:real; var k:integer; var x:real);

       var r,y: real;

begin

         k:=0;

         x:=(a+b)*0.5;

         y:=f(x);

         r:=f(a);

    while (y<>0) and (b-a>eps) do

begin

         k:=k+1;

         if r*y>0 then a:=x

     else b:=x;

          x:=(a+b)*0.5;

          y:=f(x)

      end

    end;

begin

           bisect ( a, b, eps, k, x);

     writeln('приближённое значение = ' , x:12:8);

     writeln('количество  итераций = ' , k:4);

end.

 

 

program Simpson;

uses Crt;

var

         a, b, x, Delta:               real;

         i, Step, N:                   integer;

         y, y1, y2:                    real;

         Er, H:                        real;

         ch:                           char;