Лабораторные работы по вычислительной математике

пределить графически корни уравнения:

3


. Определить графически корни уравнения:

4


. Определить графически корни уравнения:

2. Определить аналитически корни уравнения:

ЛАБОРАТОРНАЯ РАБОТА №10 «Методы решения систем линейных уравнений ». Студента группы ПВ-22 Малютина Максима.

Задание. Решить систему уравнений с точностью до 0,001. а) Методом итераций
Система :

б) Методом Ньютона. Система :

uses crt; type fun=function(x:real):real;

funcs=array[1..4] of fun;

fun2=function(x,y:real):real;

function fun_x(y:real):real; begin fun_x:=-0.4-sin(y); end;

function fun_y(x:real):real; begin fun_y:=(cos(x+1))/2; end;

function f(x,y:real):real; begin f:=sin(x+y)-1.5*x-0.1 end;

function g(x,y:real):real; begin g:=x*x+y*y-1 end;

function dfx(x,y:real):real; begin dfx:=sin(x+y)-1.5 end;

function dfy(x,y:real):real; begin dfy:=sin(x+y) end;

function dgx(x,y:real):real; begin dgx:=2*x end;

function dgy(x,y:real):real; begin; dgy:=2*y end;

Procedure Iteration(funx,funy:fun;x,y,e,q:real); var xn,yn:real;

m:byte; begin e:=abs(e*(1-q)/q); xn:=x; yn:=y; m:=0; repeat

x:=xn;y:=yn;

xn:=funx(y);

yn:=funy(x);

inc(m) until (abs(xn)+abs(yn)-abs(x)-abs(y))<e; writeln('Решение : X = ',xn,'. Y= ',yn) end;

Procedure Nuton(dfx,dfy,dgx,dgy,f,g:fun2;x,y,eps:real); var d,d1,d2,xn,yn,dx1,dy1:real; begin xn:=x;yn:=y; repeat

x:=xn;y:=yn;

d:=dfx(x,y)*dgy(x,y)-dfy(x,y)*dgx(x,y);

d1:=-f(x,y)*dgy(x,y)+g(x,y)*dfy(x,y);

d2:=-g(x,y)*dfx(x,y)+f(x,y)*dgx(x,y);

dx1:=d1/d;dy1:=d2/d;

xn:=x+dx1;

yn:=y+dy1; until (abs(xn)+abs(yn)-abs(x)-abs(y))<eps; writeln('Решение : X = ',xn,' Y= ',yn) end;

var x,y,q,eps:real; begin clrscr; writeln('Введите заданную точность'); readln(eps); writeln('Введите начальные значения X, Y '); readln(x,y); writeln('Введите q '); readln(q); Iteration(fun_x,fun_y,x,y,eps,q); writeln('Введите начальные значения X, Y '); readln(x,y); Nuton(dfx,dfy,dgx,dgy,f,g,x,y,eps) end.

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

Введите заданную точность 0.001 Введите начальные значения X, Y -0.88 0.45 Введите q 0.9 Решение : X = -8.76048170584909E-0001. Y= 4.96164420593686E-0001 Количество шагов = 7 Введите начальные значения X, Y 0.58 0.8 Решение : X = 5.89956109385639E-0001 Y= 8.07435397634436E-0001 Количество шагов = 4

ЛАБОРАТОРНАЯ РАБОТА №12 «Численное интегрирование ».

Студента группы ПВ-22 Малютина Максима.

Задание. Различными способами вычислить приближенно значение определенного интеграла.

Известно, что определенный интеграл функции типа численно представляет собой площадь криволинейной трапеции ограниченной кривыми x=0, y=a, y=b и y= (Рис. 1). Есть два метода вычисления этой площади или определенного интеграла — метод трапеций (Рис. 2) и метод средних прямоугольников (Рис. 3).

Рис. 1. Криволинейная трапеция.

Рис. 2. Метод трапеций.

Рис. 3. Метод средних прямоугольников.

По методам трапеций и средних прямоугольников соответственно интеграл равен сумме площадей прямоугольных трапеций, где основание трапеции какая-либо малая величина (точность), и сумма площадей прямоугольников, где основание прямоугольника какая-либо малая величина (точность), а высота определяется по точке пересечения верхнего основания прямоугольника, которое график функции должен пересекать в середине. Соответственно получаем формулы площадей — для метода трапеций:

,

для метода средних прямоугольников:

.

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

Остановимся поподробнее на формуле Симпсона и т.н. формуле «трех восьмых».

Формула Симпсона:

Формула «трех восьмых»:

Число разбиений n должно быть кратно трем.

Экстраполяция по Ричардсону.

Пусть In1 и In2 – два приближеных значения интуграла, найденные по одной и той же формуле при n1 и n2 (n2>n1). Тогда более точное значение этого интеграла можно найти по формуле:

,где m – порядок остаточного члена (для формулы трапеций m=2, для формулы Симпсона m=4)

Соответственно этим формулам и составим алгоритм.

Листинг программы.

program Integral;

uses Crt, Dos;

function Fx(x:real):real;

begin

fx:=(1+0.9*x*x)/(1.3+sqrt(0.5*x*x+1))

{В этом месте запишите функцию, для вычисления интеграла.}

end;

Function Yi(x,h:real;i:LongInt):real;

begin

Yi:=fx(x+i*h)

end;

Function CountBar(x1,x2,h:real):real;

var xx1,xx2:real;

c:longint;

i:real;

begin

writeln('-->Метод средних прямоугольников.');

i:=0;

for c:=1 to round(abs(x2-x1)/h) do

begin

write('Итерация ',c,chr(13));

xx1:=Fx(x1+c*h);

xx2:=Fx(x1+c*h+h);

i:=i+abs(xx1+xx2)/2*h

end;

writeln('------------------------------------------------');

CountBar:=i

end;

Function CountTrap(x1,x2,h:real):real;

var xx1,xx2,xx3:real;

c:longint;

i:real;

begin

writeln('--> Метод трапеций.');

i:=0;

for c:=1 to round(abs(x2-x1)/h) do

begin

write('Итерация ',c,chr(13));

xx1:=Fx(x1+c*h);

xx2:=Fx(x1+c*h+h);

if xx2>xx1 then xx3:=xx1 else xx3:=xx2;

i:=i+abs(xx2-xx1)*h+abs(xx3)*h

end;

writeln('------------------------------------------------');

CountTrap:=i

end;

Function CountSimpson(x1,x2,h:real):real;

var i:real;

j,n:LongInt;

begin

n:=round(abs(x2-x1)/h);

writeln('-->Метод Симпсона.');

i:=fx(x1);

j:=2;

while j<=n-1 do

begin

i:=i+4*yi(x1,h,j)+2*yi(x1,h,j+1);

j:=j+2

end;

writeln('------------------------------------------------');

CountSimpson:=h/3*(i+4*yi(x1,h,j)+yi(x1,h,j+1));

end;

Function CountThree(x1,x2,h:real):real;

var s1,s2,s3:real;

i,n:LongInt;

begin

writeln('-->Метод "Трех восьмых".');

n:=round((abs(x2-x1))/h);

if n mod 3=0 then

begin

s1:=fx(x1)+fx(x2);

s2:=0;s3:=0;

for i:=1 to n do

begin

if i mod 3=0 then s3:=s3+yi(x1,h,i)

else s2:=s2+yi(x1,h,i)

end;

CountThree:=3*h/8*(s1+3*s2+2*s3);

writeln('------------------------------------------------')

end

else writeln('Неверное число шагов !!! (Должно быть кратно 3) ')

end;

Function Richardson(i1,i2,m,a:real):double;

var b:double;

begin

b:=a/(exp(m*ln(a))-1);

Richardson:=i2+b*(i2-i1)

end;

var i1,i2,i,x1,x2,h1,h2:real;

c:byte;

n1,n2,m:word;

begin

writeln('------------------------------------------------');

writeln('-= Программа вычисления определенного интеграла =-');

writeln('Введите исходные значения: ');

write('Начальное значение x (x нижн)=');Readln(x1);

write('Конечное значение x (x верхн)=');Readln(x2);

repeat

write('Вычисление по числу итераций(1) или по шагу(2)? ');readln(c);

until (c=1) or (c=2);

case c of

1: begin

write('Количество итераций (n1)=');Readln(n1);

write('Количество итераций (n2)=');Readln(n2);

h1:=(abs(x2-x1))/n1;

h2:=(abs(x2-x1))/n2;

writeln('Шаг вычисления (h1)=',h1);

writeln('Шаг вычисления (h2)=',h2)

end;

2: begin

write('Шаг вычисления (h1)=');Readln(h1);

write('Шаг вычисления (h2)=');Readln(h2);

writeln('Количество итераций (n1)=',round(abs(x2-x1)/h1));

writeln('Количество итераций (n2)=',round(abs(x2-x1)/h2))

end;

end;

i1:=CountTrap(x1,x2,h1);

writeln('Интеграл=',i1);

i2:=CountTrap(x1,x2,h2);

writeln('Интеграл=',i2);

writeln('Экстраполирование Ричардсона для случая трапеций: ');

writeln('Интеграл = ',Richardson(i1,i2,2,n2/n1));

readln;

i1:=CountBar(x1,x2,h1);

writeln('Интеграл = ',i1);

i2:=CountBar(x1,x2,h2);

writeln('Интеграл = ',i2);

writeln('Экстраполирование Ричардсона для случая прямоугольников ');

writeln('Интеграл = ',Richardson(i1,i2,3,n2/n1));

writeln('------------------------------------------------');

i1:=CountSimpson(x1,x2,h1);

writeln('Интеграл = ',i1);

i2:=CountSimpson(x1,x2,h2);

writeln('Интеграл = ',i2);

writeln('Экстраполирование Ричардсона для случая Симпсона ');

writeln('Интеграл = ',Richardson(i1,i2,3,n2/n1));

i1:=CountThree(x1,x2,h1);

writeln('Интеграл = ',i1);

i2:=CountThree(x1,x2,h2);

writeln('Интеграл = ',i2);

writeln('Спасибо за использование программы ;-) ');

readln

end.

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

------------------------------------------------

-= Программа вычисления определенного интеграла =-

Введите исходные значения:

Начальное значение x (x нижн)=0.9

Конечное значение x (x верхн)=2.34

Вычисление по числу итераций(1) или по шагу(2)? 1

Количество итераций (n1)=4

Количество итераций (n2)=5

Шаг вычисления (h1)= 3.60000000000127E-0001

Шаг вычисления (h2)= 2.88000000000011E-0001

--> Метод трапеций.

------------------------------------------------

Интеграл= 3.21492525852918E-0003

--> Метод трапеций.

------------------------------------------------

Интеграл= 4.61840165326066E-0003

Экстраполирование Ричардсона для случая трапеций:

Интеграл = 7.73723808599729E-0003

------------------------------------------------

Интеграл = 2.53128978246764E-0003

Экстраполирование Ричардсона для случая прямоугольников

Интеграл = 3.65111028007424E-0003

------------------------------------------------

-->Метод Симпсона.

------------------------------------------------

Интеграл = 1.07491181758519E-0002

-->Метод Симпсона.

------------------------------------------------

Интеграл = 9.02681082661161E-0003

Экстраполирование Ричардсона для случая Симпсона

Интеграл = 6.76804708990304E-0003

------------------------------------------------

-->Метод "Трех восьмых".

Неверное число шагов !!! (Должно быть кратно 3)

Интеграл = 0.00000000000000E+0000

------------------------------------------------

-->Метод "Трех восьмых".

Неверное число шагов !!! (Должно быть кратно 3)

Интеграл = 0.00000000000000E+0000

------------------------------------------------

-->Метод Гаусса.

Интеграл = 1.40977850823276E-0002

------------------------------------------------

-->Метод Гаусса.

Интеграл = 1.40649829885291E-0002

Спасибо за использование программы ;-)

Задание 1. Отделить корни уравнения графически и уточнить один из них методом хорд с точностью до 0,001.

x=0.672275594. Количество шагов – 5.

Задание 2. Отделить корни уравнения аналитически и уточнить один из них методом хорд с точностью до 0,001.

x=-0.3219021. Количество шагов – 5.

То же самое методом хорд:

1. x=0.67235827. Количество шагов – 3.

2. x=-0.3222222. Количество шагов – 3.

Задание 1. Комбинированным методом хорд и касательных решить уравнение 3-ей степени, вычислив корни с точностью до 0,001.

X1=-0.810246. Количество шагов – 2.

X2= 1.168039. Количество шагов – 2.

X3=2.641798. Количество шагов – 2.

{определение корня методом хорд}

uses crt;

function fun(x:real):real;{заданная функция}

Begin

fun:=x+ln(x)/ln(10)-0.5;

End;

function fun2(x:real):real;{вторая производная заданной функции}

Begin

fun2:=-1/ln(10)/x/x;

End;

var a,b,e,e1,min,max,x,x1,n,f,f1:real;

m:byte;

BEGIN

clrscr;

writeln('Введите промежуток где возможен корень');

write('a=');readln(a);

write('b=');readln(b);

write('Введите точность E=');readln(e);

writeln('Введите m и M');

write('m=');readln(min);

write('M=');readln(max);

if fun(a)*fun2(a)>0 then

begin

n:=a;

x:=b;

end

else

begin

n:=b;

x:=a;

end;

f:=fun(n);

e1:=(e*min)/(max-min);

m:=0;

repeat

x1:=x;

f1:=fun(x1);

x:=x1-(f1*(n-x1))/(f-f1);

m:=m+1;

until e1>=abs(x-x1);

writeln('Корень =',x);

writeln(m);

END.

{определение корня методом Ньютона}

uses crt;

function fun(x:real):real;{заданная функция}

Begin

fun:=x*x*x+3*x+1;

End;

function fun1(x:real):real;{первая производная}

Begin

fun1:=3*(x*x+1);

End;

function fun2(x:real):real;{вторая производная}

Begin

fun2:=6*x;

End;

var a,b,e,e1,min,max,x,x1,n:real;

m:byte;

BEGIN

clrscr;

writeln('Введите промежуток где возможен корень');

write('a=');readln(a);

write('b=');readln(b);

write('Введите точность E=');readln(e);

writeln('Введите m и M');

write('m=');readln(min);

write('M=');readln(max);

if fun(a)*fun2(a)>0 then

begin

n:=b;

x:=a;

end

else

begin

n:=a;

x:=b;

end;

e1:=sqrt((2*min*e)/max);

m:=0;

repeat

x1:=x;

x:=x1-fun(x1)/fun1(x1);

m:=m+1;

until e1>=abs(x-x1);

writeln('Корень =',x);

writeln(m);

END.

{определение корня комбинированным методом}

uses crt;

function fun(x:real):real;{заданная функция}

Begin

fun:=x*x*x-3*x*x+2.5;

End;

function fun1(x:real):real;{первая производная}

Begin

fun1:=3*x*(x-2);

End;

function fun2(x:real):real;{вторая производная}

Begin

fun2:=6*x-6;

End;

procedure chord(n,x1:real;var x:real);{метод хорд}

Begin

x:=x1-(fun(x1)*(n-x1))/(fun(n)-fun(x1));

End;

procedure nuton(x1:real;var x:real);{метод Ньютона}

Begin

x:=x1-fun(x1)/fun1(x1);

End;

var x,a,b,e,xx,x1,xn,n,n1:real;

m:byte;

BEGIN

clrscr;

writeln('Введите промежуток где возможен корень');

write('a=');readln(a);

write('b=');readln(b);

write('Введите точность E=');readln(e);

if fun(a)*fun2(a)>0 then

begin

n:=a;x:=b;

n1:=b;x1:=a;

end

else

begin

n:=b;x:=a;

n1:=a;x1:=b;

end;

m:=0;

repeat

nuton(x1,xn);

chord(xn,x,xx);

x:=xx;

x1:=xn;

m:=m+1;

until abs(xx-xn)<=e;

writeln('Корень =',(xx+xn)/2);

writeln(m);

END.

Задание 1. Отделить корни уравнения графически и уточнить один из них методом итераций с точностью до 0,001.

X=0,213310688. Количество шагов – 3.

З

адание 2. Отделить корни уравнения аналитически и уточнить один из них методом итераций с точностью до 0,001.

X=-1,1246907. Количество шагов – 4.

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

uses crt;

function fun(x:real):real;

begin

fun:=x*x*x-0.1*x*x+0.4*x+2;

end;

function fun1(x:real):real;

begin

fun1:=3*x*x-0.2*x+0.4;

end;

var u,x,xn,q:real;

min,max:real;

a,b,e:real;

m:byte;

begin

clrscr;

writeln('Введите промежуток где возможен корень');

write('a=');readln(a);

write('b=');readln(b);

write('Введите точность E=');readln(e);

writeln('Введите m и M');

write('m=');readln(min);

write('M=');readln(max);

u:=2/(min+max);

q:=(max-min)/(max+min);

e:=abs(e*(1-q)/q);

x:=a;

m:=0;

repeat

xn:=x;

x:=xn-u*fun(xn);

m:=m+1;

until abs(x-xn)<e;

writeln('Корень =',x);

writeln(m);

end.

ЛАБОРАТОРНАЯ РАБОТА №3 «Алгебра матриц». Студента группы ПВ-22 Малютина Максима.

Задание. Обратить матрицу методом разбиения ее на произведение двух треугольных матриц.

Вариант 8.

При разбиении матрицы А на две треугольные, используются следующие формулы:

M=1..n.

Получены следующие результаты: Матрица T: Матрица R:

Матрица A-1:

program lab_3; { Лабораторная работа по вычмслительной математике N 3 Нахождение матрицы, обратной данной }

const Sizem = 10; { максимальная размерность матрицы }

type mattype = array [1..Sizem,1..Sizem] of double;

{ процедура для вывода матрицы на экран } procedure WriteMat (var m : mattype;n,n1 : byte); var k,i : byte; begin

writeln;

for k := 1 to n do

begin

for i := 1 to n1 do

write(m[k,i] : 7:3,' ');

writeln

end; end;

{ процедура ввода значений элементов матрицы } procedure inputmat (var m : mattype; var n : byte); var k,i : byte; begin

writeln;

write ('Размер матрицы = ');

readln(n);

for k := 1 to n do

for i := 1 to n do

read (m[k,i]); end;

{ процедура транспонирования матрицы } procedure Transpose (var m : mattype;n : byte); var k,i : byte;

ttt : double; begin

for k := 1 to n do

for i := k+1 to n do

begin

ttt := m[k,i];

m[k,i] := m[i,k];

m[i,k] := ttt;

end; end;

{ процедура умножения двух матриц (a*b=c) } procedure MulMat (a : mattype; ma,na : byte;

b : mattype; mb,nb : byte;

var c : mattype; var mc,nc : byte); var k,i,j : byte;

s : double; begin

if na = nb then

begin

mc := ma;

nc := nb;

for k := 1 to mc do

for j := 1 to nc do

begin

s := 0;

for i := 1 to nc do

s := s+a[k,i]*b[i,j];

c[k,j] := s;

end;

end

else

begin

writeln('Неправильный размер матрицы !');

halt

end; end;

{ процедура получения двух треугольных матриц произведение которых равно матрице m } procedure GetRnT (var m,r,t : mattype; n : byte); var k,i,m1,l : byte; begin

for k := 1 to n do

for i := 1 to n do

begin

if k=i then r[k,i] := 1

else r[k,i] := 0;

t[k,i] := 0;

end;

for m1 := 1 to n do

begin

if m1 = 1 then

begin

for i := 1 to n do

t[i,1] := m[i,1];

for i := 2 to n do

r[1,i] := m[1,i]/t[1,1];

end

else

begin

k := m1;

for i := m1 to n do

begin

t[i,k] := m[i,k];

for l := 1 to k-1 do

t[i,k] := t[i,k] - t[i,l]*r[l,k];

end;

i := m1;

for k := i+1 to n do

begin

r[i,k] := m[i,k];

for l := 1 to i-1 do

r[i,k] := r[i,k] - t[i,l]*r[l,k];

r[i,k] := r[i,k] / t[i,i];

end;

end;

end; end;

{ процедура обращения нижней треугольной матрицы } procedure BackMat (var t : mattype; n : byte); var i,k,l : byte;

x : mattype; begin

for i := 1 to n do

x[i,i] := 1/t[i,i];

for k := 1 to n-1 do

for i := k+1 to n do

begin

x[i,k] := 0;

for l := k to i-1 do

x[i,k] := x[i,k] - t[i,l]*x[l,k];

x[i,k] := x[i,k]/t[i,i];

end;

t := x end;

var m,m1,r,t : mattype;

n : byte; { ------------- основная программа ---------------- } begin

writeln ('Лабораторная работа N 2 ');

InputMat(m,n); { ввод матрицы m }

GetRnT(m,r,t,n);{получение треугольных матриц t и r}

Writeln('Матрица T: ');

WriteMat(t,n,n);

readln;

Writeln('Матрица R: ');

WriteMat(r,n,n);

readln;

BackMat(t,n); { обращение матрицы t }

Transpose(r,n); { транспонирование матрицы r }

BackMat(r,n); {обращение матрицы r (транcпонир.)}

Transpose(r,n);{транспонирование обращенной м-цы r }

MulMat(r,n,n,t,n,n,m1,n,n);

{получение матрицы,обратной матрице m}

WriteMat (m1,n,n);{ печать обратной матрицы }

readln;

MulMat(m,n,n,m1,n,n,m,n,n); { Проверка вычислений }

WriteMat(m,n,n);

readln; end.

ЛАБОРАТОРНАЯ РАБОТА №4 «Методы решения систем линейных уравнений ». Студента группы ПВ-22 Малютина Максима.

Задание. Решить систему по схеме Халецкого с точностью до 0,0001.

Вариант 8.

При разбиении матрицы А на две треугольные, используются следующие формулы:

M=1..n.

Получены следующие результаты: Матрица T: Матрица R:

Матрица X:

program lab_3; { Лабораторная работа по вычмслительной математике N 3 Нахождение матрицы, обратной данной }

const Sizem = 10; { максимальная размерность матрицы }

type mattype = array [1..Sizem,1..Sizem] of double;

{ процедура для вывода матрицы на экран } procedure WriteMat (var m : mattype;n,n1 : byte); var k,i : byte; begin

writeln;

for k := 1 to n do

begin

for i := 1 to n1 do

write(m[k,i] : 7:3,' ');

writeln

end; end;

{ процедура ввода значений элементов матрицы } procedure inputmat (var m : mattype; var n : byte); var k,i : byte; begin

writeln;

write ('Размер матрицы = ');

readln(n);

for k := 1 to n do

for i := 1 to n do

read (m[k,i]); end;

{ процедура транспонирования матрицы } procedure Transpose (var m : mattype;n : byte); var k,i : byte;

ttt : double; begin

for k := 1 to n do

for i := k+1 to n do

begin

ttt := m[k,i];

m[k,i] := m[i,k];

m[i,k] := ttt;

end; end;

{ процедура умножения двух матриц (a*b=c) } procedure MulMat (a : mattype; ma,na : byte;

b : mattype; mb,nb : byte;

var c : mattype; var mc,nc : byte); var k,i,j : byte;

s : double; begin

if na = nb then

begin

mc := ma;

nc := nb;

for k := 1 to mc do

for j := 1 to nc do

begin

s := 0;

for i := 1 to nc do

s := s+a[k,i]*b[i,j];

c[k,j] := s;

end;

end

else

begin

writeln('Неправильный размер матрицы !');

halt

end; end;

{ процедура получения двух треугольных матриц произведение которых равно матрице m } procedure GetRnT (var m,r,t : mattype; n : byte); var k,i,m1,l : byte; begin

for k := 1 to n do

for i := 1 to n do

begin

if k=i then r[k,i] := 1

else r[k,i] := 0;

t[k,i] := 0;

end;

for m1 := 1 to n do

begin

if m1 = 1 then

begin

for i := 1 to n do

t[i,1] := m[i,1];

for i := 2 to n do

r[1,i] := m[1,i]/t[1,1];

end

else

begin

k := m1;

for i := m1 to n do

begin

t[i,k] := m[i,k];

for l := 1 to k-1 do

t[i,k] := t[i,k] - t[i,l]*r[l,k];

end;

i := m1;

for k := i+1 to n do

begin

r[i,k] := m[i,k];

for l := 1 to i-1 do

r[i,k] := r[i,k] - t[i,l]*r[l,k];

r[i,k] := r[i,k] / t[i,i];

end;

end;

end; end;

{ процедура обращения нижней треугольной матрицы } procedure BackMat (var t : mattype; n : byte); var i,k,l : byte;

x : mattype; begin

for i := 1 to n do

x[i,i] := 1/t[i,i];

for k := 1 to n-1 do

for i := k+1 to n do

begin

x[i,k] := 0;

for l := k to i-1 do

x[i,k] := x[i,k] - t[i,l]*x[l,k];

x[i,k] := x[i,k]/t[i,i];

end;

t := x end;

var m,m1,r,t : mattype;

n : byte; { ------------- основная программа ---------------- } begin

writeln ('Лабораторная работа N 2 ');

InputMat(m,n); { ввод матрицы m }

GetRnT(m,r,t,n);{получение треугольных матриц t и r}

Writeln('Матрица T: ');

WriteMat(t,n,n);

readln;

Writeln('Матрица R: ');

WriteMat(r,n,n);

readln;

BackMat(t,n); { обращение матрицы t }

Transpose(r,n); { транспонирование матрицы r }

BackMat(r,n); {обращение матрицы r (транcпонир.)}

Transpose(r,n);{транспонирование обращенной м-цы r }

MulMat(r,n,n,t,n,n,m1,n,n);

{получение матрицы,обратной матрице m}

WriteMat (m1,n,n);{ печать обратной матрицы }

readln;

MulMat(m,n,n,m1,n,n,m,n,n); { Проверка вычислений }

WriteMat(m,n,n);

readln; end.

ЛАБОРАТОРНАЯ РАБОТА №5 «Методы решения систем линейных уравнений ».

Студента группы ПВ-22 Малютина Максима.

Задание. Решить систему линейных уравнений методом квадратных корней с точностью до 0,001.

Вариант 8.

При разбиении матрицы А на треугольную используются следующая формулы:

j=1..n.

const size=10;

type vector=array[1..size] of real;

matrix=array[1..size] of vector;

Procedure InputVector(var a:vector;n:byte);

var i:byte;

begin

for i:=1 to n do

begin

writeln('Введите ',i,'-ый элемент ');

readln(a[i]);

end;

end;

Procedure InputMatrix(var a:matrix;n:byte);

var i:byte;

begin

for i:=1 to n do

begin

writeln('Введите ',i,'-ую строку матрицы ');

InputVector(a[i],n)

end;

end;

Procedure OutputVector(var a:vector;n:byte);

var i:byte;

begin

for i:=1 to n do write(a[i]:10:5);

writeln

end;

Procedure OutputMatrix(var a:matrix;n:byte);

var i:byte;

begin

for i:=1 to n do outputvector(a[i],n)

end;

Procedure GetT(var t:matrix;a:matrix;n:byte);

var i,j,l:byte;

s:real;

begin

for i:=1 to n do

for j:=1 to n do t[i,j]:=0;

for j:=1 to n do

begin

s:=0;

for l:=1 to j-1 do s:=s+sqr(t[j,l]);

s:=a[j,j]-s;

t[j,j]:=sqrt(s);

for i:=j+1 to n do

begin

s:=0;

for l:=1 to j-1 do s:=s+t[i,l]*t[j,l];

t[i,j]:=(a[i,j]-s)/t[j,j]

end;

end;

end;

procedure MulMatrix(a:matrix;ma,na:byte;b:matrix;mb,nb:byte;var c:matrix;var mc,nc:byte);

var i,j,k:byte;

s:real;

begin

if na=nb then

begin

mc:=ma;

nc:=nb;

for k:=1 to mc do

for j:=1 to nc do

begin

s:=0;

for i:=1 to nc do

s:=s+a[k,i]*b[i,j];

c[k,j]:=s

end;

end

else

begin

writeln('Неверные размеры матриц !!! ');

halt

end;

end;

procedure MulVector(a:matrix;ma,na:byte;b:vector;nb:byte;var c:vector;var nc:byte);

var i,j:byte;

s:real;

begin

if na=nb then

begin

nc:=nb;

for i:=1 to nc do

begin

s:=0;

for j:=1 to nc do s:=s+a[i,j]*b[j];

c[i]:=s;

end;

end

else

begin

writeln('Неверные размеры матриц !!! ');

halt

end;

end;

Procedure TransposeMatrix(var a:matrix;n:byte);

var i,j:byte;

s:real;

begin

for i:=1 to n do

for j:=1 to n do

begin

s:=a[i,j];

a[i,j]:=a[j,i];

a[j,i]:=s

end;

end;

procedure GetY(t:matrix;b:vector;var y:vector;n:byte);

var i,k:byte;

s:real;

begin

for i:=1 to n do

begin

s:=0;

for k:=1 to i-1 do s:=s+t[i,k]*y[k];

y[i]:=(b[i]-s)/t[i,i];

end;

end;

procedure GetX(t:matrix;y:vector;var x:vector;n:byte);

var j,k:byte;

s:real;

begin

for j:=n downto 1 do

begin

s:=0;

for k:=j+1 to n do s:=s+t[k,j]*x[k];

x[j]:=(y[j]-s)/t[j,j];

end;

end;

var a,at,at2,t:matrix;

b,b2,y,x:vector;

n:byte;

begin

writeln('Введите размерность матрицы коэффициентов ');readln(n);

writeln('Введите элементы матрицы коэффициентов ');

InputMatrix(a,n);

writeln('Введите вектор свободных членов ');

InputVector(b,n);

at:=a;

TransposeMatrix(at,n);

MulMatrix(a,n,n,at,n,n,at2,n,n);

MulVector(at,n,n,b,n,b2,n);

Writeln('Пребразованная матрица А: ');

at:=at2;

outputmatrix(at,n);

Writeln('Преобразованный вектор B: ');

b:=b2;

outputvector(b,n);

writeln;

GetT(t,at,n);

Writeln('Пребразованная матрица T: ');

outputmatrix(t,n);

GetY(t,b,y,n);

writeln('Вектор Y');

outputvector(y,n);

GetX(t,y,x,n);

writeln('Вектор X');

outputvector(x,n)

end.

Пребразованная матрица А:Преобразованный вектор B:

4.97540 1.82880 1.26010 -0.14480 4.23870 -4.67000

1.82880 3.64830 -1.77800

1.26010 -1.77800 3.78260

Пребразованная матрицаT:Вектор Y

2.23056 0.00000 0.00000 -0.06492 2.48788 -1.05155

0.81988 1.72514 0.00000 Вектор X

0.56493 -1.29913 1.33256 -0.14090 0.84788 -0.78912

ЛАБОРАТОРНАЯ РАБОТА №6 «Методы решения систем линейных уравнений ».

Студента группы ПВ-22 Малютина Максима.

Задание. Решить систему линейных уравнений методом Гаусса с выбором максимального элемента по столбцу с точностью до 0,001.

Вариант 8.

При решении системы уравнений методом Гаусса используются следующие формулы:

Шаг № I: (i:=1, n-1)

Среди элементов i столбца (начиная с i-ой строки до n-ой) выбираем max по модулю элемент. Если их несколько, выбираем первый. Меняем местами i-ое уравнение и отмеченное.

Далее проводим i-ый шаг метода Гаусса:

j:=i+1,n mj = aji / aii; Вычисляем mj

Далее исключаем xi:

Вычитаем из строк i+1..n i-ую строку, помноженную на m:

k:=i+1,n

j:=1,n akj = akj - aij * mk

bk = bk – bi * mk

Д

алее осуществляется обратный ход метода Гаусса:

program gauss_max;

const size=10;

type vector=array[1..size] of real;

matrix=array[1..size] of vector;

Procedure InputVector(var a:vector;n:byte);

var i:byte;

begin

for i:=1 to n do

begin

writeln('Введите ',i,'-ый элемент ');

readln(a[i]);

end;

end;

Procedure InputMatrix(var a:matrix;n:byte);

var i:byte;

begin

for i:=1 to n do

begin

writeln('Введите ',i,'-ую строку матрицы ');

InputVector(a[i],n)

end;

end;

Procedure OutputVector(var a:vector;n:byte);

var i:byte;

begin

for i:=1 to n do write(a[i]:10:5);

writeln

end;

Procedure OutputMatrix(var a:matrix;n:byte);

var i:byte;

begin

for i:=1 to n do outputvector(a[i],n)

end;

Procedure MulVector(a:matrix;ma,na:byte;b:vector;nb:byte;var c:vector;var nc:byte);

var i,j:byte;

s:real;

begin

if na=nb then

begin

nc:=nb;

for i:=1 to nc do

begin

s:=0;

for j:=1 to nc do s:=s+a[i,j]*b[j];

c[i]:=s

end;

end

else

begin

writeln('Неверные размеры матриц !!! ');

halt

end;

end;

Procedure SwapVector(var a,b:vector);

var n:vector;

begin

n:=a;

a:=b;

b:=n

end;

Procedure Swap(var a,b:real);

var n:real;

begin

n:=a;

a:=b;

b:=n

end;

Procedure GetMaxEl(a:matrix;n,i:byte;var l:byte);

var k:byte;

max:real;

begin

max:=abs(a[i,i]);l:=i;

for k:=i to n do

if abs(a[k,i])>max then

begin

max:=abs(a[k,i]);

l:=k

end;

end;

Procedure GetAm(var a:matrix;var b:vector;n:byte);

var i,j,k,l:byte;

m:vector;

begin

for i:=1 to n-1 do

begin

GetMaxEl(a,n,i,l);

SwapVector(a[i],a[l]);

Swap(b[i],b[l]);

for j:=i+1 to n do m[j]:=a[j,i]/a[i,i];

for k:=i+1 to n do

begin

for j:=1 to n do a[k,j]:=a[k,j]-a[i,j]*m[k];

b[k]:=b[k]-b[i]*m[k]

end;

end;

end;

Procedure GetX(a:matrix;b:vector;n:byte;var x:vector);

var k,l:byte;

s:real;

begin

x[n]:=b[n]/a[n,n];

for k:=n-1 downto 1 do

begin

s:=0;

for l:=k+1 to n do s:=s+a[k,l]*x[l];

x[k]:=(b[k]-s)/a[k,k]

end;

end;

var a,am:matrix;

b,x,x2:vector;

n:byte;

begin

writeln('Введите размерность матрицы коэффициентов ');readln(n);

writeln('Введите элементы матрицы коэффициентов ');

InputMatrix(a,n);

writeln('Введите вектор свободных членов ');

InputVector(b,n);

am:=a;

GetAm(am,b,n);

writeln('Матрица Am ');

outputmatrix(am,n);

GetX(am,b,n,x);

writeln('Вектор X ');

outputvector(x,n);

MulVector(a,n,n,x,n,x2,n);

writeln('Проверка: Вектор X2 - умножение матрицы Am на X ');

outputvector(x2,n)

end.

Матрица А:Вектор B:

10.00000 6.00000 2.00000 0.00000 25.00000 8.00000 2.50000 1.50000

0.00000 6.00000 -2.00000 2.00000

0.00000 3.20000 0.40000 -1.00000

0.00000 -2.00000 -3.00000 4.00000

Матрица Am

10.00000 6.00000 2.00000 0.00000

0.00000 6.00000 -2.00000 2.00000

0.00000 0.00000 -3.66667 4.66667

0.00000 -0.00000 -0.00000 -0.20000

Вектор X

2.00000 1.00000 -0.50000 0.50000

Проверка: Вектор X2 - умножение матрицы Am на X

25.00000 8.00000 2.50000 1.50000

ЛАБОРАТОРНАЯ РАБОТА №7 «Методы решения систем линейных уравнений ».

Студента группы ПВ-22 Малютина Максима.

Задание. Составить программу, отладить ее на тестовом примере, рассмотренном на лекции.

Система :

При решении примера на лекции:

x1 = 0.526; x2 =0.628; x3 = 0.64; x4 = 1.2.

Векторы a, b, c, d.

a = {0; 2; 2; 3}

b = {5; 4.6; 3.6; 4.4}

c = {-1; -1; -0.8; 0}

d = {2; 3.3; 2.6; 7.2}

Прямой ход прогонки заключается в нахождении прогоночных коэффициентов:

Обратный ход метода прогонки заключается в нахождении неизвестных xn, xn-1, ... x1.

Он начинается с равенства: xn=bn+1;

const max=10;

type matrix=array[1..max] of real;

matrix_2=array[0..max] of real;

procedure input_matr(var a:matrix;n:byte;c:char);

var i:byte;

begin

for i:=1 to n do

begin

writeln('Введите ',i ,'-ый элемент массива ',c);

readln(a[i])

end

end;

procedure process(a,b,c,d:matrix;var x:matrix;n:byte);

var alfa,betta:matrix_2;

gamma,fi:matrix;

i:byte;

begin

betta[0]:=0;

alfa[0]:=0;

for i:=1 to n do

begin

gamma[i]:=b[i]+a[i]*alfa[i-1];

fi[i]:=d[i]-a[i]*betta[i-1];

alfa[i]:=-c[i]/gamma[i];

betta[i]:=fi[i]/gamma[i]

end;

x[n]:=betta[n];

for i:=n-1 downto 1 do x[i]:=alfa[i]*x[i+1]+betta[i]

end;

procedure out_matr_x(a:matrix;n:byte);

var i:byte;

begin

for i:=1 to n do writeln(i ,' корень уравнения равен ',a[i]:5:3)

end;

var i:byte;

a,b,c,d,x,gamma,fi:matrix;

alfa,betta:matrix_2;

n:byte;

begin

writeln('Введите размерность системы ');

readln(n);

if (n>=2) and (n<=10) then

begin

input_matr(a,n,'a');

input_matr(b,n,'b');

input_matr(c,n,'c');

input_matr(d,n,'d');

process(a,b,c,d,x,n);

out_matr_x(x,n)

end

else writeln('1< Размерность <=10 !!! ')

end.

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

1 корень уравнения равен 0.526

2 корень уравнения равен 0.628

3 корень уравнения равен 0.640

4 корень уравнения равен 1.200

ЛАБОРАТОРНАЯ РАБОТА №9 «Методы решения систем линейных уравнений ».

Студента группы ПВ-22 Малютина Максима.

Задание. Методом Зейделя решить систему линейных уравнений с точностью до 0,001.

Система :

Для решения системы уравнений методом Зейделя необходимо выполнения условия диагонального преобладания, после приведения к данному виду система имеет вид:

В
оспользуемся разложением матрицы А на В и С вида:

Далее найдем решение приближенное решение уравнения следующим способом: Правило остановки:

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

const size=10;

type vector=array[1..size] of real;

matrix=array[1..size] of vector;

norma=function(a:matrix;n:byte):real;

norma_v=function(a:vector;n:byte):real;

Procedure InputVector(var a:vector;n:byte);

var i:byte;

begin

for i:=1 to n do

begin

writeln('Введите ',i,'-ый элемент ');

readln(a[i]);

end;

end;

Procedure InputMatrix(var a:matrix;n:byte);

var i:byte;

begin

for i:=1 to n do

begin

writeln('Введите ',i,'-ую строку матрицы ');

InputVector(a[i],n)

end;

end;

Procedure OutputVector(var a:vector;n:byte);

var i:byte;

begin

for i:=1 to n do write(a[i]:10:5);

writeln

end;

Procedure OutputMatrix(var a:matrix;n:byte);

var i:byte;

begin

for i:=1 to n do outputvector(a[i],n)

end;

Procedure GetB(var b:matrix;a:matrix;n:byte);

var i,j:byte;

s:real;

begin

for i:=1 to n do

for j:=1 to n do

if i<>j then b[i,j]:=-a[i,j]/a[i,i]

else b[i,j]:=0;

end;

Procedure GetC(var c:vector;h:vector;n:byte;a:matrix);

var i:byte;

begin

for i:=1 to n do c[i]:=h[i]/a[i,i]

end;

Function Norma_1v(a:vector;n:byte):real;

var i:byte;

s:real;

begin

s:=a[1];

for i:=2 to n do if abs(a[i])>s then s:=abs(a[i]);

norma_1v:=s

end;

Function Norma_8v(a:vector;n:byte):real;

var i:byte;

s:real;

begin

s:=0;

for i:=1 to n do s:=s+abs(a[i]);

norma_8v:=s

end;

Function Norma_1(a:matrix;n:byte):real;

var s,norma:real;

i,j:byte;

begin

norma:=0;

for j:=1 to n do

begin

s:=0;

for i:=1 to n do s:=s+abs(a[i,j]);

if s>norma then norma:=s

end;

norma_1:=norma

end;

Function Norma_8(a:matrix;n:byte):real;

var s,norma:real;

i,j:byte;

begin

norma:=0;

for i:=1 to n do

begin

s:=0;

for j:=1 to n do

s:=s+abs(a[i,j]);

if s>norma then norma:=s

end;

norma_8:=norma

end;

procedure MulMatrix(a:matrix;ma,na:byte;b:matrix;mb,nb:byte;var c:matrix;var mc,nc:byte);

var i,j,k:byte;

s:real;

begin

if na=nb then

begin

mc:=ma;

nc:=nb;

for k:=1 to mc do

for j:=1 to nc do

begin

s:=0;

for i:=1 to nc do

s:=s+a[k,i]*b[i,j];

c[k,j]:=s

end;

end

else

begin

writeln('Неверные размеры матриц !!! ');

halt

end;

end;

Procedure SubMatr(a:matrix;var b:matrix;n:byte);

var i,j:byte;

begin

for i:=1 to n do

for j:=1 to n do b[i,j]:=a[i,j]-b[i,j]

end;

procedure MulVector(a:matrix;ma,na:byte;b:vector;nb:byte;var c:vector;var nc:byte);

var i,j:byte;

s:real;

begin

if na=nb then

begin

nc:=nb;

for i:=1 to nc do

begin

s:=0;

for j:=1 to nc do s:=s+a[i,j]*b[j];

c[i]:=s;

end;

end

else

begin

writeln('Неверные размеры !!! ');

halt

end;

end;

procedure MulVectorZ(a:matrix;n:byte;var b:vector);

var i,j:byte;

s:real;

begin

for i:=1 to n do

begin

s:=0;

for j:=1 to n do

s:=s+a[i,j]*b[j];

b[i]:=s;

end;

end;

Procedure SubVect(a,b:vector;var c:vector;n:byte);

var i:byte;

begin

for i:=1 to n do c[i]:=b[i]-a[i]

end;

Procedure AddVect(a:vector;var b:vector;n:byte);

var i:byte;

begin

for i:=1 to n do b[i]:=b[i]+a[i]

end;

var a,b,bn:matrix;

h,c,xr,x,xn:vector;

i,n:byte;

eps:real;

nor:norma;

norv:norma_v;

begin

writeln('Введите размерность матрицы коэффициентов ');readln(n);

writeln('Введите элементы матрицы коэффициентов ');

InputMatrix(a,n);

writeln('Введите вектор свободных членов H ');

InputVector(h,n);

writeln('Введите заданныю точность ');

readln(eps);

GetB(b,a,n);

GetC(c,h,n,a);

writeln('Матрица B: ');

OutputMatrix(b,n);

writeln('Вектор C: ');

OutputVector(c,n);

readln;

if (norma_1(b,n)<=norma_8(b,n)) and (norma_1(b,n)<>0) then

begin

nor:=norma_1;

norv:=norma_1v

end

else

begin

nor:=norma_8;

norv:=norma_8v

end;

eps:=eps*(1-nor(b,n))/nor(b,n);

for i:=1 to n do x[i]:=1;

MulVectorZ(b,n,x);

AddVect(c,x,n);

xn:=x;

MulVectorZ(b,n,xn);

AddVect(c,xn,n);

subvect(x,xn,xr,n);

while norv(xr,n)>eps do

begin

x:=xn;

MulVectorZ(b,n,xn);

AddVect(c,xn,n);

subvect(x,xn,xr,n)

end;

writeln('Значения X ');

OutputVector(x,n);

MulVector(a,n,n,x,n,c,n);

writeln('Проверка ');

OutputVector(c,n);

end.

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

Матрица B:

0.00000 0.06250 -0.11458

-0.34375 0.00000 -0.26563

-0.45946 -0.32432 0.00000

Вектор C:

-0.08333 1.26563 0.25676

Значения X

0.01836 1.30590 -0.17513

Проверка

-0.79990 8.10045 1.90065