Методы вычислительной математики
Задание 1
Решить заданную систему линейных уравнений методом простых итераций.
;
Перепишем в виде:
А=; B=; X=.
Ax=b (1) преобразуем равносильным образом к виду
X=Bx+c (2)
1 и 2 имеют одно множество решений.
X(0);
X(1)=Bx(0)+c;
X(2)=Bx(1)+c;
……………..
X(n+1)=Bx(n)+c; — Итерационная последовательность (в скобках номер итерации);
=x* — этот вектор (Х*) и будет являться решением системы 2, а значит и системы 1.
=;
Достаточное условие сходимости, метода итерации.
Если норма ||B||<1, то итерационная последовательность будет сходиться к истинному решению со скоростью геометрической прогрессии, знаменатель которой q=||B||, причем условие будет выполняться при любом начальном приближении. Чаще всего метод простой итерации производят в форме Яковлева.
Используя программу находим решение системы линейных уравнений
Задание №2.
Отделить корни уравнения, используя графико –аналитический метод. Найти корни уравнения методами бисекций, Ньютона или простых итераций. Выполнить проверку правильности найденных решений, вычислив невязки.
Задача №01.
Х4+5Х-3=0
Решение. Для локализации
Х4 = 3 – 5Х
Построив графики функций, определяем, что у решаемого уравнения имеется два корня, которые находятся в интервале -2 -1,5 и 0,5<x<1
1) -2 -1,5
Найдем корни уравнения
методом Ньютона. Для
Для функции f(x)= x4 + 5x -3 имеем:
Первая производная: f’(x)= 4x3 + 5;
Вторая производная: f’’(x)= 12x2;
f’(-2)=-32-5=-27;
f’’(-1,5)=-13.5-5=-18.5;
Производные не меняют знак внутри интервалов областей определения.
В качестве начального приближения можно выбрать правую границу интервала, для которого выполняется неравенство:
f-2)f’’(-2)0
M1=min|f’(x)|=18,5;
M2=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,
Действия продолжаются до тех пор пока 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,
Действия продолжаются до тех пор пока f(с) не станет равным 0, или разница между a и b не станет меньше заданной погрешности.
Зададим точность вычисления 0.0001
Вычисление с помощью
- Метод бисекций. Ответ: 0,73371887. Кол-во итераций -14.
Задание №3.
Используя обобщенные формулы трапеций и Симпсона вычислить определенные интегралы с заданной точностью. Проверку достижения требуемой точности проводить по правилу Рунге..
Задача №01. (Метод Трапеций)
dx
f(x) - непрерывная функция на отрезке
x0 = a -нижний предел интегрирования.
xn = b – верхний предел интегрирования.
X0=0;
Xn=1;
Предположим что Y=h(+f1+f2+…+fn-1) –верное решение интеграла.
Разобьем отрезок на 5 частей h=(1-0)/5=0,2 (n=5-кол-во разбиений)
f(x)=
f0=f(X0)=f(0)=0;
f1=f(X1)=f(0,2)=0,039479;
f2=f(X2)=f(0,4)=0,152203;
f3=f(X3)=f(0,6)=0,324252;
f4=f(X4)=f(0,8)=0,539193;
f5=f(X5)=f(1)=0,785398;
Y=0.2(+0,039479+0,152203+0,
Задача №01. (Метод Симпсона)
dx
Метод Симпсона-
график подынтегральной
Предположим Y=h/3(f0+f2m+42i-1+22i) верное решение интеграла.
Разобьем отрезок на 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,
Проверка Аналитическим методом
= -)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(x0)=y0.
Расчетные формулы:
Метод Эйлера
Метод Рунге-Кутта
Yi+1=Yi+(k1+2k2+2k3+k4);
k1=f(Xi;Yi);
k2=f(Xi+;Yi+);
k3=f(Xi+;Yi+);
k4=f(Xi+h; Yi+hk3).
Метод Эйлера-Коши более прост и требует меньшего количества вычислений, чем метод Рунге-Кутта, однако он менее точен; что является основным его недостатком.
Используя программные средства решаем уравнение.
Задача №02.
Найти приближенное решение задачи Коши y’’=f(x;y;y’), y(x0)=y0, y’(x0)=y0’или y’’’=f(x;y;y’;y’’), y(x0)=y0; y’(x0)=y0’; y’’(x0)=y0’’ методами Эйлера и Рунге-Кутта 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[
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]+
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)
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]+
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;
function f(x: real): real;
begin
f :=cos(x*x*x); {задание подынтегральной функции}
end;
function Simp (a, b: real; N:integer): real;
var y,h, x: real;
i: integer;
begin
h:=(b-a)/N;
y:=0;
x:=A;
for i:=1 to ((N-2 div 2 )+1) do
begin
y:=y+f(x)+4*f(x+h)+f(x+2*h);
x:=x+2*h;
end;
Simp:=y*h/3;
end;
begin
ClrScr;
writeln ('Приближенное вычисление определенных интегралов');
writeln ('----------------------------
writeln (' Метод Симпсона ');
writeln (' (с табулированием) ');
writeln;
write (' Введите четное число интервалов разбиений N: ');
readln (N);
write (' Введите нижний предел интегрирования a: ');
readln (a);
write (' Введите верхний предел интегрирования b: ');
readln (b);
write (' Введите шаг Delta: ');
readln (Delta);
writeln (' результат: ');
writeln (' ------------------------------
writeln ('! X ! Y ! Er ! ');
writeln (' ------------------------------
Step := Round ((b-a)/Delta)+1;
for i:=1 to Step do
begin
x :=a+Delta*(i-1);
y1:=Simp(a,x,N);
y2:=Simp(a,x,2*N);
Er:=(y2-y1)/15;
y :=y2+Er;
writeln (x:12:5, y:12:5, Er:12:5);
end;
writeln ('
Нажмите любую клавишу...
ch:=readkey;
end.
program slau;
uses crt;
var
a11,a12,a13,a14:real;
a21,a22,a23,a24:real;
a31,a32,a33,a34:real;
a41,a42,a43,a44:real;
b1,b2,b3,b4:real;
x10,x20,x30,x40:real;
x1,x2,x3,x4:real;
n1,n2,n3,n4:real;
eps:real; e:real;
iter,MaxIt:integer;
ch:char;
begin
Clrscr;
writeln('Решение системы четырх линеных уравнений');
writeln('---------------------
writeln('Метод Гауса-Зейделя')
write('Введите коэффициенты первого уравнения a11 a12 a13 a14:');
readln(a11,a12,a13,a14);
write('Введите коэффициенты второго уравнения a21 a22 a23 a24:');
readln(a21,a22,a23,a24);
write('Введите коэффициенты третьего уравнения a31 a32 a33 a34:');
readln(a31,a32,a33,a34);
write('Введите коэффициенты четвертого уравнения a41 a42 a43 a44:');
readln(a41,a42,a43,a44);
write('Введите свободные члены b1 b2 b3 b4:');
readln(b1,b2,b3,b4);
write('Введите начальное приближение x10 x20 x30 x40:');
readln(x10,x20,x30,x40);
write('Введите погрешность:');
readln(eps);
write('Введите максимальньное число итераций');
readln(MaxIt);
iter :=1;
repeat
x1 := (b1 - a12*x20 - a13*x30 - a14*x40)/a11;
x2 := (b2 - a21*x1 - a23*x30 - a24*x40)/a22;
x3 := (b3 - a31*x1 - a32*x2 - a34*x40)/a33;
x4 := (b4 - a41*x1 - a42*x2 - a43*x3)/a44;
E := sqrt (sqr(x1-x10)+sqr(x2-x20)+sqr(
x10 :=x1;
x20 :=x2;
x30 :=x3;
x40 :=x4;
iter :=iter+1;
until (E < eps) or (iter > MaxIt);
{Подсчет невязок}
n1:=abs(a11*x1+a12*x2+a13*x3+
n2:=abs(a21*x1+a22*x2+a23*x3+
n3:=abs(a31*x1+a32*x2+a33*x3+
n4:=abs(a41*x1+a42*x2+a43*x3+
writeln('Невязки');
writeln('n1= ',n1:8:5,' n2= ',n2:8:5,' n3= ',n3:8:5,' n4= ',n4:8:5);
writeln('Решения');
writeln('Решение системы четырх линеных уравнений');
writeln('X1= ',x10:5:3,' X2= ',x2:6:3,' X3= ',x3:6:3,' X4= ',x4:6:3);
writeln('количество итераций: ',iter - 1);
writeln('Нажмите любую
ch:=readkey;
end.
