Численные методы - вариант 5

Задание1. Разработать программу, которая:
вычисляет точное (аналитическое) и приближённое (численное) значения определённого интеграла;
выводит на экран график подынтегральной функции и подсчитанные значения определённого интеграла.
2. Разработать программу, которая:
находит численное решение системы дифференциальных уравнений;
создаёт файл, содержащий значения времени и фазовых координат для каждого шага интегрирования;
выводит на экран график изменения выходной величины во времени на основании описанных в файле данных.
1. Определённый интеграл
abfxdx.Таблица 1
Параметры определённого интеграла
№ варианта fxabnМетод интегрирования
5 1+lnx1 3 40 прямоугольников
n – количество интервалов разбиения диапазона интегрирования
2. Система обыкновенных дифференциальных уравнений
Таблица 2
Параметры системы уравнений
№ варианта Система уравнений Нелинейность, f∆kbk1k2T1T2ξ2TмодhМетод
5 б 2 2 ±0,05 1,5 4 0,5 0,04 - 6 2·10-3 3
Таблица 3
Система уравнений
б)
x1=1T1∙k1∙f∆-x1;x2=1T2∙k2∙x1-x2;x3=x2;y=x3,∆=k1t-x3Таблица 4
Вид нелинейности
2

Таблица 5
Метод численного решения
3
Метод Р-К 2-го порядка χ=1
ВведениеЛюбая область науки нуждается в обработке большого количества данных, что не обходится без трудоёмких математических расчетов. Многие серьёзные научные задачи требуют построения различных графиков, чтобы наглядно представить наблюдаемые процессы.
Написание программ на универсальных языках программирования, таких, как C++, Pascal ABC, Basic и др., позволило значительно облегчить сложные вычисления.
Данная работа посвящена решению задач численного интегрирования и решения систем дифференциальных уравнений, а также построения графиков функций в системе Pascal ABC.

1. Численное нахождение интеграла1.1. Листинг программыProgram Prog1_5;
// Программа вычисления определённого интеграла
uses graphABC, //Подключаем графический модуль
crt; // модуль работы с экраном в текстовом режиме
type
mas_real=array[1..100000] of real;
mas_int=array[1..100000] of integer;
const W = 800; // Размер графического окна по горизонтали
H = 600; //Размер графического окна по вертикали
// Подынтегральная функция
function f(x:real):real;
begin
f:=ln(x)+1;
end;
// Первообразная функция (аналитическое точное решение)
function p(x:real):real;
begin
p:=x*ln(x);
end;
// Процедура вычисления интеграла методом средних прямоугольников
procedure integral(a,b:real; n:integer; var s:real);
var h,x:real;
i:integer;
begin
h:=(b-a)/n;
x:=a+h/2;
s:=0;
for i:=1 to n do
begin
s:=s+f(x);
x:=x+h;
end;
s:=s*h;
end;
// Процедура расчёта табличных значений функции
// а также минимального и максимального значений по оси Y
procedure tabl(n:integer; a,b:real; var max,min:real; var x,y:mas_real);
var
i:integer;
hx:real;
begin
x[1]:=a;
hx:=(b-a)/(n-1);
for i:=1 to n do
begin
x[i]:=a+(i-1)*hx; // расчёт значений независимой переменной
y[i]:=f(x[i]); // расчёт значений зависимой переменной
end;
max:=y[1]; // первоначальное значение максимума
min:=y[1]; // первоначальное значение минимума
for i:=1 to n do begin
if y[i]>max then max:=y[i];
if y[i]<min then min:=y[i];
end;
end;
// Процедура построения графика функции
procedure grafik(n:integer; a,b,max,min:real; x,y:mas_real);
var
x0,y0,xLeft,yLeft,xRight,yRight,xn,yn,nx,ny,i:integer;
mx,my,dx,dy,num:real;
u,v:mas_int;
s:string;
begin SetWindowSize(W,H); //Устанавливаем размеры графического окна
// Координаты левой верхней границы системы координат:
xLeft:=50;
yLeft:=50;
// Координаты правой нижней границы системы координат:
xRight:=W-50;
yRight:=H-50;
// интервал по Х; a и b должно нацело делится на dx:
nx:=10; // количество засечек по оси x
dx:=(b-a)/(nx-1); // шаг подписей по оси x
ny:=10; // количество засечек по оси y
dy:=(max-min)/(ny-1); // шаг подписей по оси y
// Устанавливаем масштаб:
mx:=(xRight-xLeft)/(b-a); //масштаб по Х
my:=(yRight-yLeft)/(max-min); //масштаб по Y
// начало координат:
x0:=trunc(abs(a)*mx)+xLeft;
y0:=yRight-trunc(abs(min)*my);
// Рисуем оси координат:
line(xLeft,y0,xRight+10,y0); //ось ОХ
line(x0,yLeft-10,x0,yRight); //ось ОY
SetFontSize(12); //Размер шрифта
SetFontColor(clBlue); //Цвет шрифта
TextOut(xRight+20,y0-15,'X'); //Подписываем ось OX
TextOut(x0-10,yLeft-30,'Y'); //Подписываем ось OY
SetFontSize(8); //Размер шрифта
SetFontColor(clRed); //Цвет шрифта
{ Засечки по оси OX: }
for i:=1 to nx do
begin num:=a+(i-1)*dx; //Координата на оси ОХ
xn:=xLeft+trunc(mx*(num-a)); //Координата num в окне
Line(xn,y0-3,xn,y0+3); //рисуем засечки на оси OX
str(Num:2:2,s);
if abs(num)>1E-15 then // Исключаем 0 на оси OX
TextOut(xn-TextWidth(s) div 2,y0+10,s)
end;
{ Засечки на оси OY: }
for i:=1 to ny do
begin
num:=min+(i-1)*dy; //Координата на оси ОY
yn:=yRight-trunc(my*(num-min));
Line(x0-3,yn,x0+3,yn); //рисуем засечки на оси Oy str(num:2:2,s);
if abs(num)>1E-15 then //Исключаем 0 на оси OY
TextOut(x0-40,yn-TextHeight(s) div 2,s)
end;
TextOut(x0-10,y0+10,'0'); //Нулевая точка
{ График функции строим по точкам: }
for i:=1 to n do
begin
u[i]:=x0+round(x[i]*mx); //Координата Х в графическом окне
v[i]:=y0-round(y[i]*my); //Координата Y в графическом окне
end;
for i:=1 to n-1 do
begin
line(u[i],v[i],u[i+1],v[i+1]);
circle(u[i],v[i],3);
end;
circle(u[n],v[n],3); // отображение крайней правой точки графика
end;
// Основная программа
var n:integer;
a,b,min,max,s:real;
x,y:mas_real;
begina:=1; // левая граница интервала интегрирования
b:=3; // правая граница интервала интегрирования
n:=40; // число разбиений интервала интегрирования
tabl(n,a,b,max,min,x,y); // расчёт данных для построения графика
grafik(n,a,b,max,min,x,y); // построение графика функции
readln; // задержка до нажатия ввода - переключение в текстовый режим
ClearWindow; // Очистка графического окна
integral(a,b,n,s); // Вычисление интеграла численным методом
SetFontSize(10); //Размер шрифта
writeln('Интеграл по методу прямоугольников Ipryam = ',s:4:4);
writeln();
writeln('Точное значение интеграла Itoch = ',p(b)-p(a):4:4);
readln;
end.
1.2. Блок-схема алгоритмаВ программе использовано две функции и три процедуры.
1. Функция f(x) вычисляет значение подынтегральной функции определённого интеграла в зависимости от аргумента x и возвращает его.
2. Функция p(x) вычисляет значение первообразной от подынтегральной функции определённого интеграла в зависимости от аргумента x и возвращает его.
3. Процедура integral(a,b,n,s) в зависимости от введенных данных вычисляет значение определённого интеграла численным методом – методом средних прямоугольников и возвращает это значение. Данные для расчёта значения определённого интеграла: a – левая граница интервала интегрирования; b – правая граница интервала интегрирования; n – число точек разбиения интервала интегрирования. Данные выходные: s – возвращаемое значение определённого интеграла.
4. Процедура tabl(n,a,b,max,min,x,y) в зависимости от введённых параметров производит вычисление значений массивов независимой переменной и зависимой величины подынтегральной функции, а также минимальное и максимальное значение подынтегральной функции на рассматриваемом интервале. Данные для расчёта: n – число точек разбиения интервала интегрирования; a – левая граница интервала интегрирования; b – правая граница интервала интегрирования. Данные выходные: max – максимальное значение функции; min – минимальное значение функции; x – массив значений по оси абсцисс; y – массив значений по оси ординат.
5. Процедура grafik(n,a,b,max,min,x,y) в зависимости от введенных данных строит график с осями координат подынтегральной функции по массиву данных. Данные для построения графика: n – число точек разбиения интервала интегрирования; a – левая граница интервала интегрирования; b – правая граница интервала интегрирования; max – максимальное значение функции; min – минимальное значение функции; x – массив значений по оси абсцисс; y – массив значений по оси ординат.

Рис. 1.2.1 – Блок – схема основной программы

Рис. 1.2.2 – Блок – схема функции f(x)

Рис. 1.2.3 – Блок – схема функции p(x)

Рис. 1.2.4 – Блок – схема процедуры integral(a,b,n,s)

Рис. 1.2.5 – Блок – схема процедуры tabl(n,a,b,max,min,x,y)

Рис. 1.2.6 – Блок – схема процедуры grafik(n,a,b,max,min,x,y)

Рис. 1.2.7 – Блок – схема процедуры grafik(n,a,b,max,min,x,y)

Рис. 1.2.8 – Блок – схема процедуры grafik(n,a,b,max,min,x,y)
1.3. Результат работы программыПосле запуска созданной программы Prog1_5.pas на выполнение в графическом окне появляется график подынтегральной функции (рис. 1.3.1). После нажатия на клавишу ввода также появляются результаты решения определённого интеграла численным методом прямоугольников и аналитическое (точное) решение определённого интеграла (рис. 1.3.2).

Рис. 1.3.1 – График подынтегральной функции

Рис. 1.3.2 – Результаты решения определённого интеграла
Как видно из результатов решения (рис. 1.3.2) при заданных условиях значение определённого интеграла, вычисленного с помощью метода средних прямоугольников совпадает с точным значением до 3-го знака (при округлении) после запятой. При требованиях более точного решения определённого интеграла необходимо увеличить количество точек разбиения интервала интегрирования.
2. Численное решение системы дифференциальных уравнений2.1. Листинг программыProgram Prog2_5;
// Программа решения системы дифференциальных уравнений
uses graphABC; // Подключаем графический модуль
type matr_real=array[1..4,1..10000] of real;
matr_int=array[1..4,1..10000] of integer;
mas_real=array[1..10000] of real;
mas_int=array[1..10000] of integer;
const W=800; // Размер графического окна по горизонтали
H=600; //Размер графического окна по вертикали
// Параметры для системы ОДУ
// Значения коэффициентов системы ОДУ
k=2;
b=0.05;
k1=1.5;
k2=4;
T1=0.5;
T2=0.04;
Tmod=6; // правая граница интервала интегрирования
ht=0.002; // шаг интегрирования
// Функция дельта
function d(x:real):real;
begin
d:=k*(1-x);
end;
// Функция нелинейности
function fn(x:real):real;
begin
if x<=-b then fn:=-1;
if x>=b then fn:=1;
if (-b<x) and (x<b) then fn:=0;
end;
// Функция первого уравнения системы ОДУ
function f1(x1,x3:real):real;
begin
f1:=1/T1*(k1*fn(d(x3))-x1);
end;
// Функция второго уравнения системы ОДУ
function f2(x1,x2:real):real;
begin
f2:=1/T2*(k2*x1-x2);
end;
// Функция третьего уравнения системы ОДУ
function f3(x:real):real;
begin f3:=x;
end;
// Процедура записи данных в текстовый файл
procedure write_matrix(n:integer; t:mas_real; x:matr_real);
var
r:text;
i:integer;
begin
assign(r,'rez.txt'); // связывание файловой переменной с именем файла на диске (в текущем каталоге
rewrite(r); // создание и открытие файла на запись
for i:=1 to n do
writeln(r,t[i]:4:4,' ',x[1,i]:4:4,' ',x[2,i]:4:4,' ',x[3,i]:4:4,' ',x[4,i]:4:4);
close(r); // закрытие файла, сохранение всех ещё незаписанных данных на диск
end;
// Процедура чтения данных из текстового файла
procedure read_matrix(var nt:integer; var tt:mas_real; var xt:matr_real);
var
r:text;
i,j:integer;
begin
assign(r,'rez.txt');
reset(r);
nt:=0;
while not eof(r) do
begin
nt:=nt+1;
read(r,tt[nt],xt[1,nt],xt[2,nt],xt[3,nt],xt[4,nt]);
end;
close(r);
end;
// Процедура расчёта табличных значений функции
// а также минимального и максимального значений по оси Y
procedure tabl(var n:integer; var t:mas_real; var x:matr_real);
var
k11,k21,k31,k12,k22,k32:real;
i,j:integer;
begin
n:=trunc(Tmod/ht+1); // количество узловых точек для расчёта
t[1]:=0; // начальное значение отрезка интегрирования
x[1,1]:=0; // начальное условие для первой переменной
x[2,1]:=0; // начальное условие для второй переменной
x[3,1]:=0; // начальное условие для третьей переменной
x[4,1]:=x[3,1]; // начальное значение выходной переменной
for i:=2 to n do // расчёт системы ОДУ методом Р-К 2-го порядка с X=1
begin
t[i]:=t[i-1]+ht;
k11:=ht*f1(x[1,i-1],x[3,i-1]);
k21:=ht*f2(x[1,i-1],x[2,i-1]);
k31:=ht*f3(x[2,i-1]);
k12:=ht*f1(x[1,i-1]+k11/2,x[3,i-1]+k31/2);
k22:=ht*f2(x[1,i-1]+k11/2,x[2,i-1]+k21/2);
k32:=ht*f3(x[2,i-1]+k31/2);
x[1,i]:=x[1,i-1]+k12;
x[2,i]:=x[2,i-1]+k22;
x[3,i]:=x[3,i-1]+k32;
x[4,i]:=x[3,i];
end;
end;
// Процедура построения графика функции
procedure grafik(n:integer; t:mas_real; x:matr_real);
var
x0,y0,xk,yk,kvox,kvoy,hu,hv,nx,ny,i:integer;
a1,b1,ck,dk,gk,hk,hx,hy,min,max:real;
u,v:mas_int;
s,stroka:string;
begin SetWindowSize(W,H); //Устанавливаем размеры графического окна
// Координаты левой верхней границы системы координат:
x0:=50;
y0:=50;
xk:=50;
yk:=50;
kvox:=11;
kvoy:=11;
a1:=0;
b1:=Tmod;
{ Изображение рамки фокруг графика.}
Rectangle(x0,y0,W-xk,H-yk);
max:=x[4,1]; // первоначальное значение максимума
min:=x[4,1]; // первоначальное значение минимума
// расчёт минимума и максимума для нижнего и верхнего значений оси Y
for i:=1 to n do
begin
if x[4,i]>max then max:=x[4,i];
if x[4,i]<min then min:=x[4,i];
end;
{ Вычисление коэффициентов пересчета.}
ck:=(W-x0-xk)/(b1-a1);
dk:=x0-ck*a1;
gk:=(H-y0-yk)/(min-max);
hk:=y0-gk*max;
{ Формирование значений в экранной системе координат.}
for i:=1 to n do
begin
u[i]:=trunc(ck*t[i]+dk);
v[i]:=trunc(gk*x[4,i]+hk);
end;
{ Выбор цвета изображения графика.}
SetPenColor(clTeal);
{ Изображение графика. }
for i:=1 to n do
circle(u[i],v[i],2);
{ Определение расстояния между линиями сетки по оси OX.}
{ в "экранной" системе координат.}
hu:=trunc((W-x0-xk)/(kvox-1));
{ Определение расстояния между линиями сетки по оси OY.}
{ в "экранной" системе координат.}
hv:=trunc((H-y0-yk)/(kvoy-1));
SetPenStyle(2); // стиль линий сетки
SetPenColor(clNavy); // цвет линий сетки
{ Изображение линий сетки, параллельных оси OY.}
for i:=1 to kvox-2 do
line(x0+i*hu,y0,x0+i*hu,H-yk);
{ Изображение линий сетки, параллельных оси OX.}
for i:=1 to kvoy-2 do
line(x0,y0+i*hv,W-xk,y0+i*hv);
// Определение расстояния между линиями сетки по оси OX и OY.}
{ в "бумажной" системе координат. }
hx:=(b1-a1)/(kvox-1);
hy:=(max-min)/(kvoy-1);
SetFontColor(clBlack); //Цвет шрифта
TextOut(W-30,H-60,'X'); //Подписываем ось OX
TextOut(x0-5,y0-30,'Y'); //Подписываем ось OY
SetFontColor(clBlue); //Цвет шрифта
for i:=1 to kvox do
begin
{ Вычисление очередного значения, выводимого под осью OX, и перевод его в строковое значение. }
Str(a1+(i-1)*hx:2:2,stroka);
{ Вывод строкового значения на экран под осью OX.}
TextOut(x0+(i-1)*hu-15,H-yk div 2 - 15, stroka);
end;
for i:=1 to kvoy do
begin
{ Вычисление очередного значения, выводимого левее оси OY, и перевод его в строковое значение. }
Str(max-(i-1)*hy:1:2,stroka);
{ Вывод строкового значения на экран левее оси oy}
TextOut(x0 div 2 - 10,y0-7+(i-1)*hv, stroka);
end;
end;
// Основная программа
var n,nt:integer;
t,tt:mas_real;
x,xt:matr_real;
begintabl(n,t,x); // расчёт данных в табличном виде
write_matrix(n,t,x); // запись данных в файл
read_matrix(nt,tt,xt); // чтение данных из файла
grafik(nt,tt,xt); // построение графика функции
end.
2.2. Блок-схема алгоритмаВ программе использовано пять функций и четыре процедуры.
1. Функция d(x) вычисляет значение дельта функции и возвращает его.
2. Функция fn(x) вычисляет значение функции нелинейности и возвращает его.
3. Функция f1(x1,x3) вычисляет значение функции первого уравнения системы дифференциальных уравнений и возвращает его.
4. Функция f2(x1,x2) вычисляет значение функции второго уравнения системы дифференциальных уравнений и возвращает его.
5. Функция f3(x) вычисляет значение функции третьего уравнения системы дифференциальных уравнений и возвращает его.
6. Процедура write_matrix(n,t,x) в зависимости от рассчитанных значений размера массивов, массива времени и массива решения системы дифференциальных уравнений производит запись в текстовый файл rez.txt.
7. Процедура read_matrix(nt,tt,xt) производит чтение данных из файла rez.txt и присваивает соответствующие значения массивам tt и xt.
8. Процедура tabl(n,t,x) в зависимости от установленных параметров (в виде констант в теле программы) производит вычисление решения системы дифференциальных уравнений численным методом в виде массива значений координат времени t, матрицы значений функций решения x и максимального и минимального значений по оси ординат. Данные выходные: n – число узловых точек, в которых производятся расчёты; t – массив значений по оси абсцисс; x – матрица значений по оси ординат.
9. Процедура grafik(nt,tt,xt) в зависимости от введенных данных строит график с осями координат решения системы дифференциальных уравнений по массивам данных. Данные для построения графика: nt – число узловых точек, в которых производятся расчёты; max – максимальное значение функции; min – минимальное значение функции; tt – массив значений по оси абсцисс; xt – матрица значений по оси ординат.

Рис. 2.2.1 – Блок – схема основной программы

Рис. 2.2.2 – Блок – схема функции d(x)

Рис. 2.2.3 – Блок – схема функции fn(x)

Рис. 2.2.4 – Блок – схема функции f1(x1,x3,t)

Рис. 2.2.5 – Блок – схема функции f2(x1,x2)

Рис. 2.2.6 – Блок – схема функции f3(x)

Рис. 2.2.7 – Блок – схема процедуры write_matrix(n,t,x)

Рис. 2.2.8 – Блок – схема процедуры read_matrix(nt,tt,xt)

Рис. 2.2.9 – Блок – схема процедуры tabl(n,t,x)

Рис. 2.2.10 – Блок – схема процедуры grafik(n,t,x)

Рис. 2.2.11 – Блок – схема процедуры grafik(n,t,x)

Рис. 2.2.12 – Блок – схема процедуры grafik(n,t,x)

Рис. 2.2.13 – Блок – схема процедуры grafik(n,t,x)

Рис. 2.2.14 – Блок – схема процедуры grafik(n,t,x)
2.3. Результат работы программыПосле запуска созданной программы Prog2_5.pas на выполнение в графическом окне появляется график решения системы дифференциальных уравнений (рис. 2.3.1).
Выходная функция на графике отображается кривой, снабжённой круглыми маркерами.
Также при выполнении программы в текстовый файл rez.txt в рабочей папке программы Pascal ABC записываются массивы данных численного решения системы дифференциальных уравнений.

Рис. 2.3.1 – График решения системы дифференциальных уравнений: выходная функция выделена на рисунке с помощью круглых маркеров

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

Список использованных источниковБахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. 3-е изд., перераб. и доп. М.: БИНОМ. Лаборатория знаний, 2009.
Бахвалов Н.С. Численные методы. М.: Наука, 2006.
Справочная система программы Pascal ABC.
Цветков А.С. Язык программирования Pascal. Система программирования ABC Pascal. – Спб., 2012.