Программирование решения алгебраических уравнений

2

• Оглавление
• 1. Решение системы линейных алгебраических уравнений методом Гаусса
• 2. Решение нелинейных уравнений
• 2.1 Программа решения нелинейных алгебраических уравнений методом дихотомии (половинного деления)
• 2.2 Программа решения нелинейных алгебраических уравнений методом Ньютона (методом касательных)
• 3. Численное интегрирование
• 3.1 Программа вычисляющую определенный интеграл по формуле средних прямоугольников
• 3.2 Программа вычисляющую определенный интеграл по формуле Симпсона
• Библиографический список

• 1. Решение системы линейных алгебраических уравнений методом Гаусса
• Задача: Написать программу решения системы линейных алгебраических уравнений n-го порядка А * х = B методом Гаусса. Программа должна проверять правильность найденного решения подстановкой его в систему.
• -48,2*х1 - 4,8*х2 + 25,8*х3 - 21,6*х4 + 10,8*х5 = 36,6
• 49,4*х1 + 37*х2 - 23,4*х3 + 0,2*х4 + 40,2*х5 = -3,5
• 46,8*х1 + 2,4*х2 - 21,5*х3 - 33,8*х4 + 22,4*х5 = -48,7
• -47,2*х1 - 39,9*х2 + 9,4*х3 + 22,5*х4 - 35,1*х5 = -24,4
• 39,8*х1 + 24,8*х2 + 45,6*х3 - 41,3*х4 + 24,2*х5 = 33,5
• Описание алгоритма
• В переменную n вводится порядок матрицы системы. С помощью вспомогательной процедуры Input вводятся двумерный массив A и одномерный массив b, после чего оба массива и переменная n передаются функции Gauss. В фукции Gauss для каждого k-го шага вычислений выполняется поиск максимального элемента в k-м столбце матрицы начинаяя с k-й строки. Номер строки, содержащей максимальный элемент сохраняеется в переменной l. В том случае если максимальный элемент находится не в k-й строке, строки с номерами k и l меняются местами. Если же все эти элементы равны нулю, то происходит прекращение выполнения функции Gaussc результатом false. После выбора строки выполняется преобразование матрицы по методу Гаусса.
• Метод Гаусса:
• На 1-м шаге метода среди элементов aij определяют максимальный по модулю элемент ai1j1. Первое уравнение системы и уравнение с номером i1 меняют местами. Далее исключаем неизвестное х1 из уравнений с номерами i = 2, 3, …, n. Найдем величины qi1 = ai1/a11 (i = 2, 3, …, n), называемые множителями 1-го шага. Вычтем последовательно из второго, третьего, …, n-го уравнений системы первое уравнение, умноженное соответственно на q21, q31, …, qn1. Это позволит обратить в нуль коэффициенты при x1 во всех уравнениях, кроме первого. После ислючаем неизвестное x2 из уравнений с номерами i = 3, 4, …, n. Пусть a22(1) ? 0, где a22(1) - коэффициент, называемый главным (или ведущим) элементом 2-го шага. Вычислим множители 2-го шага qi2 = ai2(1) / a22(1) (i = 3, 4, …, n) и вычтем последовательно из третьего, четвертого, …, n-го уравнения системы второе уравнение, умноженное соответственно на q32, q42, …, qm2.
• Аналогично проводятся остальные шаги. Опишем очередной k-й шаг.
• На k-м шаге метода среди коэффициентов aij(k-1) при неизвестных в уравнениях системы с номерами i = k, …, n выбирают максимальный по модулю коэффициент aikjk(k-1). Затем k-е уравнение и уравнение, содержащее найденный коэффициент, меняют местами и исключают неизвестное xjk из уравнений с номерами i = k + 1, …, n.
• На этапе обратного хода неизвестные вычисляют в следующем порядке: xjn, xjn-1, …, xj1. Из последнего уравнения системы находим xn. Подставляя найденное значение xn в предпоследнее уравнение, получим xn-1. Осуществляя обратную подстановку, далее последовательно находим xn-1, xn-2, …, x1.
• Далее полученное решение системы и помещается в массив x. Потом с помощью процедуры Proverka решение подставляется в систему уравунений и проверяется равенство массивов х = b. Полученное решение выводится на экран при помощи вспомогательной процедуры Output.
• Текст программы:
• Program Metod_Gaussa;
• Uses crt;
• Const maxn=10;
• Type
• Data=Real;
• Matrix=Array[1..maxn,1..maxn] of Data;
• Vector=Array[1..maxn] of Data;
• Procedure Input(n:integer;var a:Matrix;var b:Vector);{Procedura vvoda matrix}
• var i,j:integer;
• begin
• for i:=1 to n do
• begin
• for j:=1 to n do
• begin
• Write('Vvedite Matrix[',i,',',j,']=');
• readln (a[i,j]);
• end;
• end;
• Writeln;
• For i:=1 to n do
• begin
• Write('Vvedite Vector b',i,'=');
• read(b[i]);
• end;
• end;
• Function Gauss(n:integer;a:Matrix;b:Vector;var x:Vector): Boolean;
• var
• i,j,k,l,e,s:integer;
• q,m,t:Data;
• begin
• for k:=1 to n-1 do
• begin {Naidem stroku s max elementom v k-tom stolbce}
• l:=0;
• m:=0;
• For i:=k to n do
• begin
• If Abs(a[i,k])>m then
• begin
• m:=Abs(a[i,k]);
• l:=i;
• end;
• end;
• {Esli y vsex strok ot k do n element v k-tom stolbce nulevoi, to sistema ne
• imeet odnoznachnogo reshenia}
• If l=0 then
• begin
• Gauss:=False;
• Exit;
• end;
• {Menyaem l stroky s k-toi}
• If l<>k then
• begin
• For j:=1 to n do
• begin
• t:=a[k,j];
• a[k,j]:=a[l,j];
• a[l,j]:=t;
• end;
• t:=b[k];
• b[k]:=b[l];
• b[l]:=t;
• end;
• end;
• {Preobrazuem matrix}
• For k:=1 to n-1 do
• begin
• For i:=k+1 to n do
• begin
• q:=a[i,k]/a[k,k];
• For j:=1 to n do
• begin
• { If J=K then a[i,j]:=0
• else}
• a[i,j]:=a[i,j]-q*a[k,j];
• end;
• b[i]:=b[i]-q*b[k];
• end;
• end;
• {Vichislyaem reshenie}
• x[n]:=b[n]/a[n,n];
• For i:=n-1 downto 1 do
• begin
• t:=0;
• For j:=1 to n-i do
• begin
• t:=t+a[i,i+j]*x[i+j];
• x[i]:=(1/a[i,i])*(b[i]-t);
• end;
• end;
• Gauss:=true;
• end;
• procedure Proverka(var n,w:integer;var x,p,b:Vector;var a:Matrix);{Proverka REZULTATA}
• var i,j:integer;
• e:real;
• begin
• e:=0;
• for i:=1 to n do
• begin
• p[i]:=0;
• for j:=1 to n do
• begin
• e:=a[i,j]*x[j];
• p[i]:=p[i]+e;
• end;
• If Abs(b[i])=Round(p[i]) then w:=w+1;
• end;
• end;
• procedure Output(var n,w:integer;var x,p,b:Vector;var a:Matrix);{Vvivod REZULTATA}
• var i,j:integer;
• begin
• Writeln('Isxodnaya matrix');
• for i:=1 to n do
• begin
• for j:=1 to n do
• begin
• Write(a[i,j]:7:2,'*x',j,' ');
• if j=n then Write('=',b[i]:6:2);
• end;
• writeln;
• end;
• Writeln;
• Writeln('REZULTAT:');
• Writeln;
• for i:=1 to n do
• begin
• Writeln('x',i,'=',x[i]:5:2);
• end;
• Writeln;
• Writeln('Proverka (Podstavim REZULTAT v sistemy yranenei)');
• Writeln;
• for i:=1 to n do
• begin
• for j:=1 to n do
• begin
• Write('(',a[i,j]:5:2,'*',x[j]:5:2,')');
• if j<>n then write('+');
• end;
• Write('=',p[i]:5:2);
• Writeln;
• end;
• if w=n then writeln('Sistema reshena VERNO!')
• else writeln('Sistema reshena NE VERNO!!!');
• end;
• Var n,i,w:integer;
• a:Matrix;
• b,x,p:Vector;
• begin
• ClrScr;
• Writeln('Vvedite poryadoc matrix(max=10)');
• Repeat Read(n);
• Until (n>0) and (n<=maxn);
• ClrScr;
• Writeln('Vvedite matrix');
• Writeln;
• Input(n,a,b);
• If Gauss(n,a,b,x)=true then
• begin
• Clrscr;
• Proverka(n,w,x,p,b,a);
• Output(n,w,x,p,b,a);
• end
• Else Writeln('Dannyiu sistemy nevozmozno reshit');
• Readkey;
• end.
• Программа решения систем линейных уравнений по методу Гаусса
• Vvedite poryadoc matrix(max=10)
• 5
• Vvedite matrix
• -48,2 -4,8 25,8 -21,6 10,8
• 49,4 37 -23,4 0,2 40,2
• 46,8 2,4 -21,5 -33,8 22,4
• -47,2 -39,9 9,4 22,5 -35,1
• 39,8 24,8 45,6 -41,3 24,2
• Vvedite vektor
• 36,6
• -3,5
• 33,5
• -24,4
• -48,7
• Isxodnaya matrix
• -48.20*x1 -4.80*x2 25.80*x3 -21.60*x4 10.80*x5 = 36.60
• 49.40*x1 37.00*x2 -23.40*x3 0.20*x4 40.20*x5 = -3.50
• 46.80*x1 2.40*x2 -21.50*x3 -33.80*x4 22.40*x5 =-48.70
• -47.20*x1 -39.90*x2 9.40*x3 22.50*x4 -35.10*x5 =-24.40
• 39.80*x1 24.80*x2 45.60*x3 -41.30*x4 24.20*x5 = 33.50
• REZULTAT:
• x1=-0.83
• x2= 2.30
• x3= 0.35
• x4=-0.42
• x5=-0.98
• Proverka (Podstavim REZULTAT v sistemy yranenei)
• (-48.20*-0.83)+(-4.80* 2.30)+(25.80* 0.35)+(-21.60*-0.42)+(10.80*-0.98)=36.60
• (49.40*-0.83)+(37.00* 2.30)+(-23.40* 0.35)+( 0.20*-0.42)+(40.20*-0.98)=-3.50
• (46.80*-0.83)+( 2.40* 2.30)+(-21.50* 0.35)+(-33.80*-0.42)+(22.40*-0.98)=-48.70
• (-47.20*-0.83)+(-39.90* 2.30)+( 9.40* 0.35)+(22.50*-0.42)+(-35.10*-0.98)=-24.40
• (39.80*-0.83)+(24.80* 2.30)+(45.60* 0.35)+(-41.30*-0.42)+(24.20*-0.98)=33.50
• Sistema reshena VERNO!!!

• 2. Решение нелинейных уравнений
• 2.1 Задача: Написать программу решения нелинейных алгебраических уравнений методом дихотомии (половинного деления).
• Х² * 2^х = 1
• Описание алгоритма
• Пусть непрерывная функция f(x) на концах отрезка [a,b] имеет значения разных знаков, т.е. f(a)?f(b) < 0, тогда на отрезке имеется хотя бы один корень.
• Возьмем середину отрезка х=(a+b)/2. Если f(a)?f(х) < 0, то корень явно принадлежит отрезку от a до (a+b)/2 и в противном случае от (a+b)/2 до b. Поэтому берем подходящий из этих отрезков, вычисляем значение функции в его середине и т.д. до тех пор, пока длина очередного отрезка не окажется меньше заданной предельной точности (b-a) < е. Выводим х.
• Текст программы:
• Program Metod_polovinogo_deleniya;
• uses crt;
• var a,b,e,x:Real;
• function F(x:Real):Real;
• begin
• F:=x*x*Exp(x*ln(2))-1;
• end;
• begin
• ClrScr;
• Write('Vvedite gran a=');
• Readln(a);
• Write('Vvedite gran b=');
• Readln(b);
• Write('Vvedite e=');
• Readln(e);
• While (b-a)>=e do
• begin
• x:=(a+b)/2;
• if F(a)*F(x)<=0 then b:=x
• else a:=x;
• end;
• Writeln('Koren yravneniya raven: ',x:7:4);
• Readkey;
• end.
• Программа решения нелинейных уравнений методом дихотомии
• Vvedite gran a= 2
• Vvedite gran b= 3
• Vvedite e= 0.001
• Koren yravneniya raven: 2.9990
• 2.2 Задача: Написать программу решения нелинейных алгебраических уравнений методом Ньютона (методом касательных).
• Х³ + 3х² + 6х - 1 = 0
• Описание алгоритма
• Дана непрерывная функция f(x), которая содержит единственный корень на отрезке [a,b], где b>a при чем определены непрерывны и сохраняют знак f`(x) f``(x). Точность е.
• Выбираем грубое приближение корня х₀. Найдем значение функции точке х₀ и проведем касательную до пересечения с осью абсцисс, получим значение х₁. Определим значение функции в точке х₁, через эту точку проводим касательную и получаем точку х₂. Повторим этот процесс n раз
• x_n = x_n-1 - f(x_n-1) / f'(x_n-1), пока |f(x)| < e или |x_n - x_n-1| < e и получим приближенный корень уравнения.
• Текст программы:
• Program Newton;
• uses crt;
• var x0,xn,xn1,E:Real;
• it:integer;
• function F(x:Real):Real;
• begin
• F:=(x*x*x+3*x*x+6*x-1)/(3*x*x+6*x+6);
• end;
• begin
• ClrScr;
• Write('Vvedite nachalnoe pribligenie ');
• Readln(x0);
• Write('Vvedite tochnoct ');
• Readln(E);
• xn1:=x0;
• xn:=1e10;
• it:=0;
• While Abs(xn1-xn)>E do
• begin
• xn:=xn1;
• xn1:=xn-F(xn);
• it:=it+1;
• end;
• Writeln('Koren raven ',xn:7:4);
• Writeln('Kolichestvo iterachiy= ',it);
• Readkey;
• end.
• Программа решения нелинейных уравнений методом Ньютона (методом касательных)
• Vvedite nachalnoe pribligenie 3
• Vvedite tochnoct 0.001
• Koren raven 0.1542
• Kolichestvo iterachiy = 6

• 3. Численное интегрирование
• 3.1 Задача: Написать программу вычисляющую определенный интеграл по формуле средних прямоугольников
• Описание алгоритма
• При вычислении определенного интеграла следует помнить, каков геометрический смысл определенного интеграла. Если f(x)>=0 на отрезке [a; b], то определенный интеграл численно равен площади фигуры, ограниченной графиком функции y=f(x), отрезком оси абсцисс, прямой x=a и прямой x=b (рис. 1.1) Таким образом, вычисление интеграла равносильно вычислению площади криволинейной трапеции.
• Разделим отрезок [a; b] на n равных частей, т.е. на n элементарных отрезков. Длина каждого элементарного отрезка h. Точки деления будут: x=a; x1=a+h; x2=a+2*h, ... , xn-1=a+(n-1)*h; xn=b. Площадь криволинейной трапеции приближенно заменяется площадью многоугольника, составленного из n прямоугольников. Таким образом, вычисление определенного интеграла сводится к нахождению суммы n элементарных прямоугольников.
• Вводим нижний и верхний пределы функции. Вводим n(количество отрезков). Вычисляем длину отрезков h=(b-a)/n и находим точку х(середина отрезка). Затем последовательно вычисляем значения функции в точках х и суммируем результаты. Выводим сумму.
• Текст программы:
• program Metod_Pryamoygolnikov;
• uses crt;
• var a,b,s,xi,x,h:Real;
• n,i:integer;
• Function F(x:real):Real;
• begin
• F:=1/(Sqrt(x*x+2.5));
• end;
• begin
• ClrScr;
• Write('Vvedite gran integrala a=');
• Readln(a);
• Write('Vvedite gran integrala b=');
• Readln(b);
• Write('Vvedite n=');
• Readln(n);
• h:=(b-a)/n;
• x:=a+h/2;
• S:=0;
• For i:=1 to n-1 do
• begin
• xi:=x+h*i;
• S:=S+F(xi)*h;
• end;
• Write('Znachenie integrala ',S:7:5);
• Readkey;
• end.
• Программа вычисляющая определенный интеграл по формуле средних прямоугольников
• Vvedite gran integrala a=1.6
• Vvedite gran integrala b=2.2
• Vvedite n=20
• Znachenie integrala 0.22992
• 3.2 Задача: Написать программу вычисляющую определенный интеграл по формуле Симпсона.
• Описание алгоритма
• Вводим нижний и верхний пределы функции. Вводим n(количество отрезков). Делим отрезок [a, b] на n частей и получим длину каждого отрезка по формуле h = (b - a)/n. Затем по циклу будем находить точку xi и считать значение функции в этой точке. Далее находим площадь криволинейной трапеции по формуле Симпсона.
• Текст программы:
• program Metod_Simpson;
• uses crt;
• var a,b,s,h,xi:Real;
• n,i,mp:integer;
• Function F(x:real):Real;
• begin
• F:=(x*x+1)*sin(x-0.5);
• end;
• begin
• ClrScr;
• Write('Vvedite gran integrala a=');
• Readln(a);
• Write('Vvedite gran integrala b=');
• Readln(b);
• Write('Vvedite n=');
• Readln(n);
• h:=(b-a)/n;
• S:=0;
• mp:=1;
• For i:=1 to n-1 do
• begin
• xi:=a+i*h;
• S:=S+(3+mp)*F(xi);
• mp:=-mp;
• end;
• S:=h*(F(a)+F(b)+S);
• Write('Znachenie integrala ',S:7:5);
• Readkey;
• end.
• Программа вычисляющая определенный интеграл по формуле Симпсона
• Vvedite gran integrala a=0.8
• Vvedite gran integrala b=1.6
• Vvedite n=10
• Znachenie integrala 3.98249
• 
• Библиографический список
• 1. Курс высшей математики: Учебник для вузов / В.С. Шипачев; Под ред. А.Н. Тихонова. - 3-е изд. - М.: Оникс, 2007. - 600с.
• 2. Турбо Паскаль 7.0 Самоучитель. - СПб.: Питер; К.: Издательская группа BHV, 2002. - 416с.