Вычислительные методы

3.5.2 Метод равномерного перебора

 

 

 

 

 

 - +

 

 

 

 

 - +

 

 

 

 

 

 

 

 

 

 

 

 

3.6 Текст программы (написан в пакете Turbo Pascal 7.0)

3.6.1 Метод «золотого» сечения

Program mn;

uses crt;

const EPS=0.001;

Var a,b,x,Fa,Fb,Fx:real;

    k:integer;

Function F(x:real):real;

begin F:=ln(x)-x+1.8=0; end;

Begin

Clrscr;

Repeat

Writeln('please, inter A and B');

Readln(a,b);

Fa:=F(a);

Fb:=F(b);

if Fa*Fb>0 then writeln('please,inter other A and B');

until Fa*Fb<0;

k:=0;

While ABS(b-a)>EPS do begin

                     X:=(a+b)/2;

                     Writeln(x:8:4);

                     k:=k+1;

                     Fx:=F(x);

                     if Fx=0 then writeln('kopeHb=',x);

                     if Fa*Fx<0 then begin

                                     b:=x;

                                     Fb:=Fx;

                                     end

                     else begin

                          a:=x;

                          Fa:=Fx;

                          end;

                      end;

Writeln('KopeHb=',x:8:5);

Writeln('F(x)=',Fx);

Writeln('k=',k);

Readkey;

End.

 

 

3.6.2 Метод равномерного перебора.

 

Program EXTRE;

uses crt;

var

A,B,Xmin,Ymin,x,h,y:real;

function f(x:real):real;

begin

f:=ln(x)*ln(x)/x;

end;

begin

clrscr;

writeln('vvedite interval i shag ckanirovania');

readln(a,b,h);

Ymin:=1E20;

x:=A;

writeln('       x           y');

while x<b do begin

Y:=f(x);

if Y<Ymin then

 

begin

Ymin:=Y;

Xmin:=X;

end;

writeln(x:8:4,'   ', y:8:4);

x:=x+h;

end;

writeln(' min=',Xmin:8:4);

writeln('znachtnie min=',Ymin:8:4);

readln;

end.

 

 

 

 

 

 

 

 

 

 

3.7 Вывод на экран

3.7.1 Метод «золотого»  сечения

 

3.7.2 Метод  равномерного перебора

 

3.9 Вывод

Значение точки минимума, полученное в результате работы в пакете MathCad 14.0 совпадает со значениями, полученными в результате работы в пакете TurboPascal 7.0  с помощью метода «золотого» сечения и метода равномерного перебора. Метод «золотого» сечения при оптимизации выгоднее, так как при нем нет вероятности того, что точка минимума будет пропущена программой из-за задания слишком большого шага.

4. Задача 4 (вычисление определенного  интеграла)

4.1 Постановка задачи

Вычислить значение интеграла на заданном отрезке интегрирования [a;b]. Считать заданным число разбиений отрезка интегрирования n и численный метод решения. На печать вывести приближенное, точное значение интеграла и относительную погрешность вычисления(точность вычисления).

4.2 Исходные данные

 

4.3 Математическая модель

4.3.1 Метод Гаусса

Пусть система линейных алгебраических уравнений  имеет n неизвестных и ранг ее матрицы  равен n, то есть система имеет единственное решение.

Предположим, что а_11№ 0 и более того, что а₁₁ - максимальный коэффициент по модулю в матрице А. Этого всегда можно добиться перестановкой уравнений, а также переименованием переменных, то есть перестановкой столбцов.

Разделим  обе части первого уравнения  на а₁₁. Получим x₁+с₁₂х₂+…+с_1nх_n=d₁ (4₁)

Умножим обе части полученного уравнения  на коэффициент а₂₁ и вычтем его из второго уравнения. Слагаемое, содержащее х₁ во втором уравнении исчезнет.

Далее умножим уравнение 4 на а₃₁ и вычтем его из третьего уравнения системы, х₁ тоже пропадет. И так далее.

В результате получим систему, порядок которой  на единицу меньше порядка исходной системы, так как в уравнениях системы отсутствуют слагаемые  с х₁ и отсутствует первое уравнение вида (4₁)

 

 

 

4.3.2 Метод  трапеций

Аппроксимация в этом методе осуществляется полиномом первой степени. Суть метода ясна из рисунка. 

На единичном интервале

.

 
В случае равномерной сетки (h = const )

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

(5)

 

 

 

 

 

 

 

 

4.4 Решение задачи в  пакете MathCad 14.0

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

4.5 Блок-схема

4.5.1 Метод Гаусса

 

 

 

 

 

 

 

4.5.2 Метод трапеций

 

 -

 +

 

 

 

4.6 Текст программы (написан  в программе TurboPascal 7.0)

4.6.1 Метод Гаусса

Program INTEGRAL_GAUSS8;

uses crt;

var a,b,integ:real;

function f(x:real):real;

begin

f:=ln(x)*ln(x)*x;

end;

function gauss8 (a,b:real):real;

var a1,a2,x,g:real;

i:integer;

ag,xg:array [1..8] of real;

begin

ag[1]:=0.10122854;

ag[2]:=0.22238104;

ag[3]:=0.31370664;

ag[4]:=0.36268378;

ag[5]:=0.36268378;

ag[6]:=0.31370664;

ag[7]:=0.22238104;

ag[8]:=0.10122854;

xg[1]:=(-0.96028986);

xg[2]:=(-0.79666648);

xg[3]:=(-0.52553242);

xg[4]:=(-0.18343464);

xg[5]:=0.18343464;

xg[6]:=0.52553242;

xg[7]:=0.79666648;

xg[8]:=0.96028986;

a:=0;

b:=4;

a1:=(a+b)*0.5;

a2:=(b-a)*0.5;

g:=0.0;

for i:=1 to 8 do

begin

x:=a1-a2*xg[i];

g:=g+ag[i]*f(x);

end;

gauss8:=g*a2;

end;

begin

clrscr;

integ:=gauss8(a,b);

writeln('integ=  ',integ:8:2);

readkey;

end.

 

4.6.2 Метод трапеций

program Trap;

uses crt;

const p=3.14;

var a,b,x,dx,z,dz,zt:real;

    n,i,k:integer;

function f(x:real):real;

begin

f:=ln(x)*ln(x)/x;

end;

begin

clrscr;

k:=0;

a:=1;

b:=4;

zt:=1;

n:=60;

while k<5 do

begin

z:=(f(a)+f(b))/2;

dx:=(b-a)/n;

x:=a;

for i:=1 to n-1 do

begin

x:=x+dx;

z:=z+f(x);

end;

z:=z*dx;

writeln('  deleniy=   znachenia integrala');

writeln('    ',n,'    ',z:8:2,'    ',zt:8:2);

k:=k+1;

n:=n*2;

readkey;

end.

 

4.7 Вывод на экран

4.7.1 Метод Гаусса

4.7.2 Метод трапеций

4.8 Вывод

Значение  интеграла, полученное в результате работы в пакете MathCad 14.0 почти совпадает со значениями, полученными в результате работы в пакете TurboPascal 7.0 с помощью метода Гаусса и метода трапеций. Метод трапеций ищет значение интеграла точнее.

 

5. Задача 5 (Поиск частного решения дифференциального уравнения)

5.1 Постановка задачи

Найти частный интеграл (частное  решение) дифференциального уравнения.

5.2 Исходные данные

Дифференциальное уравнение: y’=y*e^x²;

Начальные условия: y(0)=1.

5.3 Математическая модель

5.3.1 Задача Коши

Задача Коши для дифференциального  уравнения n-го порядка – это нахождение такого решения этого уравнения, которое будет удовлетворять заданным начальным условиям.

5.3.2 Метод Эйлера

Наиболее простой численный  метод решения (систем) обыкновенных дифференциальных уравнений. Впервые  описан Леонардом Эйлером в 1768 году в работе «Интегральное исчисление»^[1]. Метод Эйлера является явным, одношаговым методом первого порядка точности, основанном на аппроксимации интегральной кривой кусочно линейной функцией, т. н. ломаной Эйлера.

Описание метода

Пусть дана задача Коши для уравнения  первого порядка

где функция f определена на некоторой области . Решение разыскивается на интервале [x₀,b). На этом интервале введем узлы

Приближенное решение в узлах x_i, которое обозначим через y_i определяется по формуле

Эти формулы обобщаются на случай систем обыкновенных дифференциальных уравнений.

Оценка погрешности

Метод Эйлера является методом первого  порядка. Если функция f непрерывна в D и непрерывно дифференцируема по переменной y в D, то имеет место следующая оценка погрешности

где h — средний шаг, то есть существует C > 0 такая, что .

Заметим, что условия гладкости  на правую часть, гарантирующие единственность решения задачи Коши, необходимы для  обоснования сходимости метода Эйлера.

Значение метода Эйлера

Метод Эйлера являлся исторически  первым методом численного решения  задачи Коши. О. Коши использовал этот метод для доказательства существования  решения задачи Коши. В виду не высокой  точности и вычислительной неустойчивости для практического нахождения решений  задачи Коши метод Эйлера применяется  редко. Однако в виду своей простоты метод Эйлера находит свое применение в теоретических исследованиях  дифференциальных уравнений, задач  вариационного исчисления и ряда других математических проблем.

Модифицированный метод Эйлера с пересчетом

Вычисления по методу Эйлера с пересчетом делаются в два этапа.

Прогноз:

.

Коррекция:

.

Модифицированный метод Эйлера с пересчетом имеет второй порядок  точности, однако для его реализации необходимо дважды вычислять правую часть функции. Заметим, что метод  Эйлера с пересчетом представляет собой  разновидность методов Рунге-Кутта (предиктор-корректор).

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

Метод позволяет  решать системы обыкновенных дифференциальных уравнений (ОДУ) первого порядка  следующего вида:

которые имеют решение:

где t - независимая  переменная (например, время); X, Y и т.д. - искомые функции (зависимые от t переменные). Функции f, g и т.д. - заданы. Также предполагаются заданными  и начальные условия, т.е. значения искомых функций в начальный  момент.

Одно диф. уравнение - частный случай системы с одним  элементом. Поэтому, далее речь пойдет для определенности о системе  уравнений.

Метод может быть полезен и для решения диф. уравнений высшего (второго и  т.д.) порядка, т.к. они могут быть представлены системой диф. уравнений  первого порядка.

Метод Рунге-Кутта  заключается в рекурентном применении следующих формул:

 

 

 

5.4 Решение задачи в  пакете MatchCad 14.0

5.4.1 Стандартный метод решения с использованием блока Odesolve

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

5.4.2 Метод Рунге-Кутта с использованием блока rkfixed

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

5.5 Блок-схема

5.5.1 Модифицированный метод Эйлера

 

 

 

 

 

 

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

 

 

 

 

 

5.6. Текст программы (написан  в пакете Turbo Pascal 7.0)

5.6.1 Модифицированный метод Эйлера

program dddd;

   uses crt;

   var

   x0,y0,h,x,y,xn:real;

   i,n:integer;

   function fd(x,y:real):real;

   begin

   fd:=y*exp(x);

   end;

   begin

   clrscr;

   x0:=0;

   y0:=1;

   xn:=1;

   h:=0.1;

   n:=trunc((xn-x0)/h);

   x:=x0;

   y:=y0;

   writeln('znachenie x  znachenie y');

   writeln(x:6:2,'   ',y:8:2);

   for i:=1 to n do

   begin

   y:=y+h*fd(x,y);

   x:=x+h;

   writeln(x:6:2,'   ',y:8:2);

   end;

   readln;

   end.

 

 

 

 

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

program sl;

uses crt;

var

x0,y0,h,x,y,xn,k0,k1,k2,k3,t:real;

i,n:integer;

function fd(x,y:real):real;

begin

fd:=y*exp(x);

end;

begin

clrscr;

x0:=0;

y0:=1;

xn:=1 ;

h:=0.1;

n:=trunc((xn-x0)/h);

x:=x0;

y:=y0;

writeln('znachenie x  znachenir y');

writeln(x:6:2,'   ',y:8:6);

for i:=1 to n do

begin

k0:=h*fd(x,y);

t:=x+h/2;

k1:=h*fd(t,y+0.5*k0);

k2:=h*fd(t,y+0.5*k1);

x:=x+h;

k3:=h*fd(x,y+k2);

y:=y+(k0+2*(k1+k2)+k3)/6;

writeln (x:6:2,'   ',y:8:6);

end;

readln;

end.

 

 

5.7 Вывод на экран

5.7.1 Модифицированный метод  Эйлера

5.7.2 Метод Рунге-Кутта 5.9 Вывод

Значения дифференциального уравнения, полученные в результате работы в  пакете MathCad 14.0 почти совпадают со значениями, полученными в результате работы в пакете TurboPascal 7.0 с помощью модифицированного метода Эйлера и метода Рунге-Кутта. Метод Эйлера применяется в основном для ориентировочных расчетов. Это схема второго порядка точности метода Рунге-Кутта, поэтому при решении подобных задач целесообразнее использовать именно метод Рунге-Кутта.