Просмтор этого файла доступен только зарегистрированным пользователям. Но у нас супер быстрая регистрация: достаточно только электронной почты!

Текст из файла (страница 2)

Из этой системы видно, что , где В - нижняя треугольная матрица с диагональными элементами , равными нулю, а С - верхняя треугольная матрица с диагональными элементами, отличными от нуля, α=B+C. Следовательно

При таком способе приведения исходной СЛАУ к эквивалентному виду метод простых итераций носит название метода Якоби.

откуда

Таким образом, метод Зейделя является методом простых итераций с матрицей правых частей α=(E-B)-1C и вектором правых частей (E-B)-1β.

Разрешим систему относительно неизвестных при ненулевых диагональных элементах , aii≠0, i= 1,n (если какой-либо коэффициент на главной диагонали равен нулю, достаточно соответствующее уравнение поменять местами с любым другим уравнением). Получим следующие выражения для компонентов вектора β и матрицы α эквивалентной системы

В качестве нулевого приближения x⁽⁰⁾ вектора неизвестных примем вектор правых частей x⁽⁰⁾=β. Тогда метод простых итераций примет вид:

Из (1.19) видно преимущество итерационных методов по сравнению, например, с рассмотренным выше методом Гаусса. В вычислительном процессе участвуют только произведения матрицы на вектор, что позволяет работать только с ненулевыми элементами матрицы, значительно упрощая процесс хранения и обработки матриц.

Имеет место следующее достаточное условие сходимости метода простых итераций.

Метод простых итераций (1.19) сходится к единственному решению СЛАУ при любом начальном приближении x⁽⁰⁾, если какая-либо норма матрицы α эквивалентной системы меньше единицы

Если используется метод Якоби (выражения (1.18) для эквивалентной СЛАУ), то

достаточным условием сходимости является диагональное преобладание матрицы A, т.е.

(для каждой строки матрицы A модули элементов, стоящих на главной диагонали, больше суммы модулей недиагональных элементов). Очевидно, что в этом случае ||α||_c меньше единицы и, следовательно, итерационный процесс (1.19) сходится.

Приведем также необходимое и достаточное условие сходимости метода простых итераций. Для сходимости итерационного процесса (1.19) необходимо и достаточно, чтобы спектр матрицы α эквивалентной системы лежал внутри круга с радиусом, равным единице.

При выполнении достаточного условия сходимости оценка погрешности решения на k- ой итерации дается выражением:

Следует подчеркнуть, что это неравенство дает завышенное число итераций k, поэтому редко используется на практике

1.4 Численный метод решения задачи

Разрешим систему относительно неизвестных при ненулевых диагональных элементах , aii≠0, i= 1,n (если какой-либо коэффициент на главной диагонали равен нулю, достаточно соответствующее уравнение поменять местами с любым другим уравнением). Получим следующие выражения для компонентов вектора β и матрицы α эквивалентной системы:

При таком способе приведения исходной СЛАУ к эквивалентному виду метод простых итераций носит название метода Якоби.

В качестве нулевого приближения x⁽⁰⁾ вектора неизвестных примем вектор правых частей x⁽⁰⁾=β. Тогда метод простых итераций примет вид:

Из (1.19) видно преимущество итерационных методов по сравнению, например, с рассмотренным выше методом Гаусса. В вычислительном процессе участвуют только произведения матрицы на вектор, что позволяет работать только с ненулевыми элементами матрицы, значительно упрощая процесс хранения и обработки матриц.

Имеет место следующее достаточное условие сходимости метода простых итераций.

Метод простых итераций (1.19) сходится к единственному решению СЛАУ при любом начальном приближении x⁽⁰⁾, если какая-либо норма матрицы α эквивалентной системы меньше единицы

Если используется метод Якоби (выражения (1.18) для эквивалентной СЛАУ), то

достаточным условием сходимости является диагональное преобладание матрицы A, т.е.

(для каждой строки матрицы A модули элементов, стоящих на главной диагонали, больше суммы модулей недиагональных элементов). Очевидно, что в этом случае ||α||_c меньше единицы и, следовательно, итерационный процесс (1.19) сходится.

Приведем также необходимое и достаточное условие сходимости метода простых итераций. Для сходимости итерационного процесса (1.19) необходимо и достаточно, чтобы спектр матрицы α эквивалентной системы лежал внутри круга с радиусом, равным единице.

При выполнении достаточного условия сходимости оценка погрешности решения на k- ой итерации дается выражением:

где x^·- точное решение СЛАУ.

Процесс итераций останавливается при выполнении условия , где εε≤)(kε - задаваемая вычислителем точность.

Принимая во внимание, что из (1.20) следует неравенство , можно получить априорную оценку необходимого для достижения заданной точности числа итераций. При использовании в качестве начального приближения вектора β такая оценка определится неравенством:

откуда получаем априорную оценку числа итераций k при ||α||<1

Следует подчеркнуть, что это неравенство дает завышенное число итераций k, поэтому редко используется на практике.

1.6 Текст программы

program Yakobi;

uses crt;

const

maxn=100;

type

matrix=array[1..maxn,1..maxn] of real;

vector=array[1..maxn] of real;

vector1=array[1..maxn] of real;

var

i,j,n,k1: integer;

e,norma:real;

a: matrix;

b: vector;

x2,x3: vector1;

imya,dannble_i_rezultat,ekran:string;

procedure input(var kolvo:integer; var pogreshnostb:real; var matr1:matrix; var matr2:vector);

{Ввод исходных данных}

var i,j,code:integer;

a_s: string;

b_s: string;

begin

writeln('введите количество уравнений');

readln(kolvo);

writeln('введите погрешность вычисления');

readln(pogreshnostb);

writeln('введите матрицу коэффициентов при неизвестных');

for i:=1 to kolvo do

for j:=1 to kolvo do

begin

repeat

begin

writeln('введите элемент [',i,',',j,'] и нажмите Enter');

readln(a_s);

if (i=j) and (a_s='0') then

repeat

begin

writeln('Элементы на главной диагонали должны быть отличны от нуля. Повторите ввод');

readln(a_s);

end;

until a_s<>'0';

val(a_s,matr1[i,j],code);

end;

until code=0;

end;

writeln('введите вектор свободных коэффициентов');

for j:=1 to kolvo do

begin

repeat

writeln('введите элемент ',j,' и нажмите Enter');

readln(b_s);

val(b_s,matr2[j],code);

until code=0;

end;

end;

1.7 Тестовый пример

2. Полиномиальная интерполяция функции методом Ньютона с разделенными разностями

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

Разработать схему алгоритма и написать программу на языке Turbo Pascal 7.0 для интерполирования функции, заданной в узлах, используя метод Ньютона с разделенными разностями.

2.2 Математическая формулировка задачи

Дана табличная функция:

i	xi	yi
0	x0	y
1	x1	0
2	x2	y
...	...	1
n	xn	y
		2
		...
		y
		n

или

Точки с координатами (x_i, y_i) называются узловыми точками или узлами.

Количество узлов в табличной функции равно N=n+1.

Необходимо найти значение этой функции в промежуточной точке, например, x=D, причем .

Для решения задачи строим интерполяционный многочлен.

2.3 Обзор существующих численных методов решения задачи

Интерполяция по Лагранжу

Интерполяционный многочлен может быть построен при помощи специальных интерполяционных формул Лагранжа, Ньютона, Стерлинга, Бесселя и др.

Интерполяционный многочлен по формуле Лагранжа имеет вид:

Докажем, что многочлен Лагранжа является интерполяционным многочленом, проходящим через все узловые точки, т.е. в узлах интерполирования x_i выполняется условие L_n(x_i) = y_i. Для этого будем последовательно подставлять значения координат узловых точек таблицы в многочлен (2.1). В результате получим:

если x=x₀, то L_n(x₀) = y₀,

если x=x₁, то L_n(x₁) = y₁,

……………

если x=x_n, то L_n(x_n) = y_n.

Это достигнуто за счет того, что в числителе каждой дроби при соответствующем значении у_j, j=0,1,2,…,n отсутствует сомножитель (x-x_i), в котором i=j, а знаменатель каждой дроби получен заменой переменной х на соответствующее значение х_j.

Таким образом, интерполяционный многочлен Лагранжа приближает заданную табличную функцию, т.е. L_n(x_i) = y_i и мы можем использовать его в качестве вспомогательной функции для решения задач интерполирования, т.е. .

Чем больше узлов интерполирования на отрезке [x₀,x_n] , тем точнее интерполяционный многочлен приближает заданную табличную функцию, т.е. тем точнее равенство:

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

При решении задачи экстраполирования функции с помощью интерполяционного многочлена вычисление значения функции за пределами отрезка [x₀,x_n] обычно производят не далее, чем на один шаг h, равный наименьшей величине

так как за пределами отрезка [x₀,x_n] погрешности, как правило, увеличиваются.

Интерполяция по Ньютону

Интерполяционный многочлен по формуле Ньютона имеет вид:

(2.2)

где n – степень многочлена,

- разделенные разности 0-го, 1-го, 2-го,…., n-го порядка, соответственно.

Сплайн-интерполяция

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

Используя теорию изгиба бруса при малых деформациях, можно показать, что сплайн – это группа кубических многочленов, в местах сопряжения которых первая и вторая производные непрерывны. Такие функции называются кубическими сплайнами. Для их построения необходимо задать коэффициенты, которые единственным образом определяют многочлен в промежутке между данными точками.

Например, для некоторых функций (рис.) необходимо задать все кубические функции q₁(x), q₂(x), …q_n(x).

В наиболее общем случае эти многочлены имеют вид:

где k_ij - коэффициенты, определяемые описанными ранее условиями, количество которых равно 4n. Для определения коэффициентов kij необходимо построить и решить систему порядка 4n.

Первые 2n условий требуют, чтобы сплайны соприкасались в заданных точках:

Следующие (2n-2) условий требуют, чтобы в местах соприкосновения сплайнов были равны первые и вторые производные:

Система алгебраических уравнений имеет решение, если число уравнений соответствует числу неизвестных. Для этого необходимо ввести еще два уравнения. Обычно используются следующие условия:

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

Полученный таким образом сплайн называется естественным кубическим сплайном. Найдя коэффициенты сплайна, используют эту кусочно-гладкую полиноминальную функцию для представления данных при интерполяции.

2.4 Численный метод решения задачи

Значения f(x₀), f(x₁), … , f(x_n) , т.е. значения табличной функции в узлах, называются разделенными разностями нулевого порядка (k=0).

Отношение называется разделенной разностью первого порядка (k=1) на участке [x₀, x₁] и равно разности разделенных разностей нулевого порядка на концах участка [x₀, x₁], разделенной на длину этого участка.

Для произвольного участка [x_i, x_i+1] разделенная разность первого порядка (k=1) равна

Отношение называется разделенной разностью второго порядка (k=2) на участке [x₀, x₂] и равно разности разделенных разностей первого порядка, разделенной на длину участка [x₀, x₂].

Для произвольного участка [x_i, x_i+2] разделенная разность второго порядка (k=2) равна

Таким образом, разделенная разность k-го порядка на участке [x_i, x_i+k] может быть определена через разделенные разности (k-1)-го порядка по рекуррентной формуле:

(2.3)

Где n – степень многочлена.

Максимальное значение k равно n. Тогда i =0 и разделенная разность n-го порядка на участке [x₀,x_n] равна

,

т.е. равна разности разделенных разностей (n-1)-го порядка, разделенной на длину участка [x₀,x_n].

Разделенные разности

являются вполне определенными числами, поэтому выражение (2.2) действительно является алгебраическим многочленом n-й степени. При этом в многочлене (2.2) все разделенные алгебраическим многочленом n-й степени. При этом в многочлене (2.2) все разделенные разности определены для участков [x₀, x_0+k], .

Лемма: алгебраический многочлен (2.2), построенный по формулам Ньютона, действительно является интерполяционным многочленом, т.е. значение многочлена в узловых точках равно значению табличной функции

Докажем это. Пусть х=х₀ , тогда многочлен (2.2) равен

Пусть х=х₁, тогда многочлен (2.2) равен

Пусть х=х₂, тогда многочлен (2.2) равен

Заметим, что решение задачи интерполяции по Ньютону имеет некоторые преимущества по сравнению с решением задачи интерполяции по Лагранжу. Каждое слагаемое интерполяционного многочлена Лагранжа зависит от всех значений табличной функции y_i, i=0,1,…n. Поэтому при изменении количества узловых точек N и степени многочлена n (n=N-1) интерполяционный многочлен Лагранжа требуется строить заново. В многочлене Ньютона при изменении количества узловых точек N и степени многочлена n требуется только добавить или отбросить соответствующее число стандартных слагаемых в формуле Ньютона (2.2). Это удобно на практике и ускоряет процесс вычислений.

2.5 Схема алгоритма

На рисунке 2.1 представлена схема алгоритма решения задачи №2.

На рисунке 2.2 представлена схема алгоритма ввода исходных данных (подпрограмма-процедура Vvod).

На рисунке 2.3 представлена схема алгоритма интерполяции функции по методу Ньютона с разделенными разностями (newt)

На рисунке 2.4 представлена схема алгоритма записи данных и результата в файл (подпрограмма-процедура zapisb_v_fail).

На рисунке 2.5 представлена схема алгоритма вывода содержимого записанного файла на экран (подпрограмма-процедура outputtoscreen).

2.6 Текст программы

program newton;

uses crt,graph;

const c=10;

type matr=array[0..c,0..c] of real;

mas=array[0..c] of real;

var x,y,koef_polinoma:mas;

a:matr;

b:mas;

d1:real;

n:integer;

fail,fail1,ekran:text;

procedure Vvod(var kolvo:integer; var uzel,fun:mas);

{Процедура осуществляет ввод данных:пользователь вводит с клавиатуры

узлы интерполяции и значения функции в них. Также определяется количество узлов.}

var code,i:integer; s:string;

begin

writeln('введите количество узлов');

readln(kolvo);

kolvo:=kolvo-1;

for i:=0 to kolvo do

begin

repeat

writeln('введите ',i,'-й узел интерполирования');

readln(s);

val(s,uzel[i],code);

until code=0;

repeat

writeln('введите значение функции, соответствующее данному узлу');

readln(s);

val(s,fun[i],code);

until code=0;

end;

end;

procedure newt(var kolvo:integer; D:real; var koef,uzel,fun:mas)

var L,P:real;

begin

L:=fun[0];

P:=1;

for i:=1 to kolvo do

begin

P:=P*(D-uzel[i-1]);

for j:=1 to kolvo-i do

begin

fun[j]:=(fun[j-1]-fun[j])/(uzel[j+i]-uzel[i])

end;

koef[i]:=fun[0];

L:=L+P*fun[0];

end;

end;

procedure newt(var kolvo:integer; D:real; var koef,uzel,fun:mas) {процедура интерполяции функции методом Ньютона}

var L,P:real;

begin

L:=fun[0];

P:=1;

for i:=1 to kolvo do

begin

P:=P*(D-uzel[i-1]);

for j:=1 to kolvo-i do

begin

fun[j]:=(fun[j-1]-fun[j])/(uzel[j+i]-uzel[i])

end;

koef[i]:=fun[0];

L:=L+P*fun[0];

end;

end;

procedure zapisb(koef:mas; uzel,fun:mas; kolvo:integer; var f:text);

{В данной процедуре осуществляется запись в файл данных и результата}

var i:integer;

begin

assign(f,'interpol.txt');

rewrite(f);

for i:=0 to kolvo do writeln(f,'x= ',uzel[i]:8:4,' f(x)=',fun[i]:8:4);

writeln(f,'Интерполяционный полином');

write(f,'p(x)=',koef[0]:8:4);

for i:=1 to kolvo do if i>1 then write (f,'+(',koef[i]:8:4,')*x^',i)

else write (f,'+(',koef[i]:8:4,')*x');

close(f);

end;

procedure vblvod(var f1,f2:text);

{Вывод содержимого записанного файла на экран}

var s1:string;

begin

clrscr;

assign(f1,'interpol.txt');

reset(f1);

assigncrt(f2);

rewrite(f2);

while not eof(f1) do

assigncrt(f2);

rewrite(f2);

while not eof(f1) do

begin

Readln(f1,s1);

writeln(f2,s1);

end;

close(f2);

close(f1);

end;

procedure grafik(kolvo:integer; uzlbl,funktsiya:mas; c:mas);

{Построение графика полученной функции}

var driver,mode,Err,a1,b1,z,i,j:integer; s:string; xt,yt:real;

begin

driver:=detect;

InitGraph(driver,mode,'d:\tp7\bp\bgi');

err:=graphresult;

if err<>grok then writeln('Ошибка при инициализации графического режима')

else

begin

Setcolor(9);

line(320,0,320,480);

line(0,240,640,240);

settextstyle(smallfont,horizdir,3);

setcolor(10);

outtextxy(320,245,'0');

a1:=0;

b1:=480;

z:=-10;

for i:=0 to 20 do

begin

if z<>0 then

begin

str(z,s);

setcolor(10);

outtextxy(a1,245,s);

outtextxy(325,b1,s);

setcolor(8);

line(0,b1,640,b1);

line(a1,0,a1,480);

end;

outtextxy(325,b1,s);

setcolor(8);

line(0,b1,640,b1);

line(a1,0,a1,480);

end;

a1:=a1+32;

b1:=b1-24;

z:=succ(z);

end;

setcolor(5);

for i:=0 to kolvo do

begin

xt:=uzlbl[i];

yt:=funktsiya[i];

putpixel(round(320+xt*32),round(240-yt*24),5);

end;

moveto(round(320+uzlbl[0]*32),round(240-funktsiya[0]*24));

setcolor(11);

for i:=0 to kolvo do

begin

xt:=uzlbl[i];

yt:=0;

for j:=0 to kolvo do yt:=yt+c[j]*vozvedenie_v_stepenb(uzlbl[i],j);

lineto(round(320+xt*32),round(240-yt*24));

moveto(round(320+xt*32),round(240-yt*24));

end;

readln;

closegraph;

end;

end;

{Основная часть программы}

BEGIN

CLRSCR;

Writeln('Программа осуществляет интерполирование функции, заданной в узлах');

Vvod(N,X,Y);

writeln(‘введите значение промежуточной точки’);

readln(d1);

2.7 Тестовый пример

writeln('Нажмите Enter');

readln;

newt(N,d1,X,Y,koef_polinoma);

zapisb(koef_polinoma,x,y,n,fail);

vblvod(fail,fail1);

writeln('Нажмите Enter для просмотра графика функции, затем еще раз для выхода из программы');

readln;

grafik(N,X,Y,koef_polinoma);

END.

readln(d1);

writeln('Нажмите Enter');

readln;

newt(N,d1,X,Y,koef_polinoma);

zapisb(koef_polinoma,x,y,n,fail);

vblvod(fail,fail1);

writeln('Нажмите Enter для просмотра графика функции, затем еще раз для выхода из программы');

readln;

grafik(N,X,Y,koef_polinoma);

END.

Дана табличная функция:

Вычислить разделенные разности 1-го, 2-го, 3-го порядков (n=3) и занести их в диагональную таблицу.

Разделенные разности первого порядка:

Разделенные разности второго порядка:

Разделенная разность третьего порядка:

Интерполяционный многочлен Ньютона для заданной табличной функции имеет вид:

График интерполяционного многочлена будет таким:

procedure zapisb(koef:mas; uzel,fun:mas; kolvo:integer; var f:text);

{В данной процедуре осуществляется запись в файл данных и результата}

var i:integer;

begin

assign(f,'interpol.txt');

rewrite(f);

for i:=0 to kolvo do writeln(f,'x= ',uzel[i]:8:4,' f(x)=',fun[i]:8:4);

writeln(f,'Интерполяционный полином');

write(f,'p(x)=',koef[0]:8:4);

for i:=1 to kolvo do if i>1 then write (f,'+(',koef[i]:8:4,')*x^',i)

else write (f,'+(',koef[i]:8:4,')*x');

close(f);

end;

Характеристики

Список файлов курсовой работы