бесплано рефераты

Разделы

рефераты   Главная
рефераты   Искусство и культура
рефераты   Кибернетика
рефераты   Метрология
рефераты   Микроэкономика
рефераты   Мировая экономика МЭО
рефераты   РЦБ ценные бумаги
рефераты   САПР
рефераты   ТГП
рефераты   Теория вероятностей
рефераты   ТММ
рефераты   Автомобиль и дорога
рефераты   Компьютерные сети
рефераты   Конституционное право
      зарубежныйх стран
рефераты   Конституционное право
      России
рефераты   Краткое содержание
      произведений
рефераты   Криминалистика и
      криминология
рефераты   Военное дело и
      гражданская оборона
рефераты   География и экономическая
      география
рефераты   Геология гидрология и
      геодезия
рефераты   Спорт и туризм
рефераты   Рефераты Физика
рефераты   Физкультура и спорт
рефераты   Философия
рефераты   Финансы
рефераты   Фотография
рефераты   Музыка
рефераты   Авиация и космонавтика
рефераты   Наука и техника
рефераты   Кулинария
рефераты   Культурология
рефераты   Краеведение и этнография
рефераты   Религия и мифология
рефераты   Медицина
рефераты   Сексология
рефераты   Информатика
      программирование
 
 
 

Численные методы решения типовых математических задач

Численные методы решения типовых математических задач

ФЕДЕРАЛЬНОЕ АГЕНСТВО ПО ОБРАЗОВАНИЮ

ГОСУДАРСТВЕННОЕ УЧРЕЖДЕНИЕ ВЫСШЕГО

ПРОФЕССИОНАЛЬНОГО ОБРАЗОВАНИЯ

Тульский государственный университет

Кафедра автоматики и телемеханики

Численные методы решения типовых математических задач

Курсовая работа

по дисциплине «Вычислительная математика»

Выполнил студент группы _____________

(подпись)

Проверил _____________

(подпись)

Тула 2008г.


Содержание

Введение

1. Решение систем линейных алгебраических уравнений методом простой итерации

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3. Среднеквадратическое приближение функции

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

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

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

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

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

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

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

4. Численное интегрирование функций методом Гаусса

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

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

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

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

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

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

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

Заключение

Список использованных источников


Введение

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

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

Разработка и исследование вычислительных алгоритмов и их применение к решению конкретных задач составляет содержание огромного раздела современной математики — вычислительной математики.

Вычислительную математику определяют в широком смысле этого термина как раздел математики, включающий круг вопросов, связанных с использованием ЭВМ, и в узком смысле — как теорию численных методов и алгоритмов решения поставленных математических задач.

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

При решении задачи па ЭВМ мы всегда получаем не точное решение исходной задачи, а некоторое приближенное решение. Чем же обусловлена возникающая погрешность? Можно выделить три основные причины возникновения погрешности при численном решении исходной математической задачи. Прежде всего, входные данные исходной задачи (начальные и граничные условия, коэффициенты и правые части уравнений) всегда задаются с некоторой погрешностью. Погрешность численного метода, обусловленную неточным заданием входных данных, принято называть неустранимой погрешностью. Далее, при замене исходной задачи дискретной задачей возникает погрешность, называемая погрешностью дискретизации или, иначе, погрешностью метода. Наконец, конечная разрядность чисел, представляемых в ЭВМ, приводит к ошибкам округления, которые могут нарастать в процессе вычислений

Численные методы дают приближенное решение задачи. Это значит, что вместо точного решения и (функции или функционала) некоторой задачи мы находим решение у другой задачи, близкое в некотором смысле (например, по норме) к искомому. Основная идея всех методов — дискретизация или аппроксимация (замена, приближение) исходной задачи другой задачей, более удобной для решения на ЭВМ, причем решение аппроксимирующей задачи зависит от некоторых параметров, управляя которыми, можно определить решение с требуемой точностью. Например, в задаче численного интегрирования такими параметрами являются узлы и веса квадратурной формулы. Далее, решение дискретной задачи является элементом конечномерного пространства.


1. Решение систем линейных алгебраических уравнений методом простой итерации

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

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

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

Пусть А – невырожденная матрица image036001и нужно решить систему

где диагональные элементы матрицы А ненулевые.

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

Метод Гаусса

В методе Гаусса матрица СЛАУ с помощью равносильных преобразований преобразуется в верхнюю треугольную матрицу, получающуюся в результате прямого хода. В обратном ходе определяются неизвестные.

Пусть дана СЛАУ

Запишем расширенную матрицу системы:


На первом шаге алгоритма Гаусса выберем диагональный элемент (если он равен 0, то первую строку переставляем с какой-либо нижележащей строкой) и объявляем его a11≠0 ведущим, а соответствующую строку и столбец, на пересечении которых он стоит - ведущими. Обнулим элементы ведущего столбца. Для этого сформируем числа (-a22/a11), (-a31/a11), .. , (an1/a11).

LU-разложения матриц

Компьютерная реализация метода Гаусса часто осуществляется с использованием LU-разложения матриц.

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

A=LU

где L - нижняя треугольная матрица (матрица, у которой все элементы, находящиеся выше главной диагонали равны нулю, lij=0 при i>j), U- верхняя треугольная матрица (матрица, у которой все элементы, находящиеся ниже главной диагонали равны нулю, uij=0 при i>j).

LU – разложение может быть построено с использованием описанного выше метода Гаусса. Рассмотрим k - ый шаг метода Гаусса, на котором осуществляется обнуление поддиагональных элементов k - го столбца матрицы. Как было описано выше, с этой целью используется следующая операция

В терминах матричных операций такая операция эквивалентна умножению A(k)=MkA(k-1), где элементы матрицы определяются следующим образом

В терминах матричных операций такая операция эквивалентна умножению A(k)=MkA(k-1), где элементы матрицы определяются следующим образом

В результате прямого хода метода Гаусса получим , A(n-1)=U

где A(n-1)=U - верхняя треугольная матрица, а - нижняя треугольная матрица, имеющая вид .

Таким образом, искомое разложение A=LU получено.

Метод прогонки

Метод прогонки является одним из эффективных методов решения СЛАУ с трех - диагональными матрицами, возникающих при конечно-разностной аппроксимации задач для обыкновенных дифференциальных уравнений (ОДУ) и уравнений в частных производных второго порядка и является частным случаем метода Гаусса. Рассмотрим следующую СЛАУ:

решение которой будем искать в виде

где Qi,Pi,i=1,n - прогоночные коэффициенты, подлежащие определению. Для их определения выразим из первого уравнения СЛАУ (1.1) x1 через x2, получим:

откуда

Из второго уравнения СЛАУ (1.1) с помощью (1.3) выразим x2 через x3, получим:

откуда

Продолжая этот процесс, получим из i-го уравнения СЛАУ (1.1):

следовательно

Из последнего уравнения СЛАУ имеем

то есть

Метод простых итераций

При большом числе уравнений прямые методы решения СЛАУ (за исключением метода прогонки) становятся труднореализуемыми на ЭВМ прежде всего из-за сложности хранения и обработки матриц большой размерности. В то же время характерной особенностью ряда часто встречающихся в прикладных задачах СЛАУ является разреженность матриц. Число ненулевых элементов таких матриц мало по сравнению с их размерностью. Для решения СЛАУ с разреженными матрицами предпочтительнее использовать итерационные методы.

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

Метод Зейделя решения СЛАУ

Метод простых итераций довольно медленно сходится. Для его ускорения существует метод Зейделя, заключающийся в том, что при вычислении компонента xik+1 вектора неизвестных на (k+1)-й итерации используются x1k+1, x2k+1, .., xik+1 уже вычисленные на (k+1)-й итерации. Значения остальных компонент берутся из предыдущей итерации. Так же как и в методе простых итераций строится эквивалентная СЛАУ и за начальное приближение принимается вектор правых частей . Тогда метод Зейделя для известного вектора на k-ой итерации имеет вид:

предыдущей итерации. Так же как и в методе простых итераций строится эквивалентная СЛАУ и за начальное приближение принимается вектор правых частей . Тогда метод Зейделя для известного вектора на k-ой итерации имеет вид:

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

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

откуда

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

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

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

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

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

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

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

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

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

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

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

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

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

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


Untitled-1 копи

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

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

Untitled-1 копи

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

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

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

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

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

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

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

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

Untitled-1 копи

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

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

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

Untitled-1 копи

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

Untitled-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;

L_n(x)=\frac{(x-x_1)(x-x_2)(x-x_3) \ldots \ldots (x-x_n)}{(x_0-x_1)(x_0-x_2)(x_0-x_3) \ldots (x_0-x_n)} \cdot y_0 +\\ + \frac{(x-x_0)(x-x_2)(x-x_3) \ldots \ldots (x-x_n)}{(x_1-x_0)(x_1-x_2)(x_1-x_3) \ldots (x_1-x_n)} \cdot y_1 +\\ + \frac{(x-x_0)(x-x_1)(x-x_3) \ldots \ldots (x-x_n)}{(x_2-x_0)(x_2-x_1)(x_2-x_3) \ldots (x_2-x_n)} \cdot y_2 + \ldots\\ + \frac{(x-x_0)(x-x_1)(x-x_1) \ldots \ldots (x-x_n-1)}{(x_n-x_0)(x_n-x_1)(x_n-x_1) \ldots (x_n-x_n-1)} \cdot y_n.

 

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

или

y_i = f(x_i), i=\overline{0,n}.Точки с координатами (xi, yi) называются узловыми точками или узлами.

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

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

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

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

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

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

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

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

если x=x0, то Ln(x0) = y0,

если x=x1, то Ln(x1) = y1,

……………

если x=xn, то Ln(xn) = yn.

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

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

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


f(x_k) \approx L_n(x_k).

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

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

\left|x_{i+1} – x_i\right|,

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

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

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

L_n(x) = f(x_0) + (x – x_0) \cdot f(x_0;x_1) +\\ + (x – x_0) \cdot(x – x_1) \cdot f(x_0;x_1;x_2) +\\ + (x – x_0) \cdot(x – x_1) \cdot(x – x_2) \cdot f(x_0;x_1;x_2;x_3) + \ldots +\\ + (x – x_0) \cdot(x – x_1) \cdot \ldots \cdot (x – x_{n-1}) \cdot f(x_0;x_1; \ldots;x_n),(2.2)

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

f(x_0), f(x_0;x_1), f(x_0;x_1;x_2), f(x_0;x_1;\ldots; x_n)- разделенные разности 0-го, 1-го, 2-го,…., n-го порядка, соответственно.


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

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

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

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

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

q_i(x) =k_{1j} + k_{2i} x + k_{3i} x^2 + k_{4i} x^3, i=\overline{1,n},

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

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

q_i(x_i) = y_i, i=\overline{1,n};\\q_{i+1}(x_i) = y_i, i=\overline{0,n-1}.

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

q_{i+1}'(x_i) = q'_i(x_i),&amp; i=\overline{1,n-1};\\q''_{i+1}(x_i) = q_i''(x_i),&amp; i=\overline{0,n-1}.

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

q_1^''(x_0) = 0; q_n^''(x_n) = 0.

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

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

 

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

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

Отношение f(x_0;x_1) = \frac{f(x_1)-f(x_0)}{x_1 – x_0}называется разделенной разностью первого порядка (k=1) на участке [x0, x1] и равно разности разделенных разностей нулевого порядка на концах участка [x0, x1], разделенной на длину этого участка.

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

f(x_i;x_{i+1}) = \frac{f(x_{i+1})-f(x_i)}{x_{i+1} – x_i}.

Отношение f(x_0;x_1;x_2) = \frac{f(x_1;x_2)-f(x_0;x_1)}{x_2 – x_0}называется разделенной разностью второго порядка (k=2) на участке [x0, x2] и равно разности разделенных разностей первого порядка, разделенной на длину участка [x0, x2].

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

f(x_i;x_{i+1};x_{i+2}) = \frac{f(x_{i+1};x_{i+2}) - f(x_i;x_{i+1})}{x_{i+2} – x_i}

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

f(x_i;x_{i+1};x_{i+2}; \ldots;x_{i+k}) = \frac{ f(x_{i+1};x_{i+2}; \ldots;x_{i+k}) -  f(x_i;x_{i+1}; \ldots;x_{i+k-1})}{x_{i+k} – x_i}(2.3)

Где k=\overline{1,n}, k=\overline{0,n-k}, n – степень многочлена.

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

f(x_0;x_i; \ldots;x_n) = \frac{ f(x_1;x_2; \ldots;x_n) - f(x_0;x_1; \ldots;x_n-1)}{x_n – x_0},

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

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

f(x_0;x_1), f(x_0;x_1;x_2), \ldots, f(x_0;x_1;\ldots; x_n)

являются вполне определенными числами, поэтому выражение (2.2) действительно является алгебраическим многочленом n-й степени. При этом в многочлене (2.2) все разделенные алгебраическим многочленом n-й степени. При этом в многочлене (2.2) все разделенные разности определены для участков [x0, x0+k], k=\overline{1,n}.

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

L_n(x_i) = f(x-i) = y_i; i=0,1, \ldots n.

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

L_n(x_0) = f(x_0) = y_0.

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

L_n(x_1) = f(x_0) + (x_1 – x_0) \cdot f(x_0;x_1) = f(x_0) + (x_1 – x_0) \cdot \frac{f(x_1) – f(x_0)}{x_1 – x_0} = f(x_1) = y_1

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

L_n(x_1) = f(x_0) + (x_2 – x_0) \cdot f(x_0;x_1) + (x_2 – x_0) \cdot (x_2 – x_1) \cdot f(x_0;x_1;x_2) =\\= f(x_0) + (x_2 – x_0) \cdot \frac{f(x_1) – f(x_0)}{x_1 – x_0}+  (x_2 – x_0) \cdot (x_2 – x_1) \cdot \frac{f(x_1;x_2) – f(x_0;x_1)}{x_2 – x_0}= f(x_1) = y_1

Заметим, что решение задачи интерполяции по Ньютону имеет некоторые преимущества по сравнению с решением задачи интерполяции по Лагранжу. Каждое слагаемое интерполяционного многочлена Лагранжа зависит от всех значений табличной функции yi, 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) и занести их в диагональную таблицу.

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

f(x_0;x_1) = \frac{f(x_1) – f(x_0)}{x_1 – x_0}= \frac{1,098613 – 0,693147}{3 - 2} = 0,405466.\\f(x_1;x_2) = \frac{f(x_2) – f(x_1)}{x_2 – x_1}= \frac{1,386295 – 1,098613}{4 - 3} = 0,287682.\\f(x_2;x_3) = \frac{f(x_3) – f(x_2)}{x_3 – x_2}= \frac{1,609438 – 1,386295}{5 - 4} = 0,223143.

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

f(x_0;x_1;x_2) = \frac{f(x_1;x_2) – f(x_0;x_1)}{x_2 – x_0}= \frac {0,287682 – 0,405466}{4 - 2} = -0,058892.\\f(x_1;x_2;x_3) = \frac{f(x_2;x_3) – f(x_0;x_2)}{x_3 – x_1}= \frac {0,223143 – 0,287682}{5 - 3} = - 0,0322695.

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

f(x_0;x_1;x_2;x_3) = \frac{f(x_1;x_3) – f(x_0;x_2)}{x_3 – x_0}= \frac {- 0,0322695 – (- 0,058892)}{5 - 2} = 0,00887416

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

L_3(x) = f(x_0) + (x – x_0) \cdot f(x_0;x_1) + (x-x_0)(x  - x_1) \cdot f(x_0;x_1;x) +\\+ (x – x_0)(x – x_1)(x – x_2) \cdot f(x_0;x_1;x_2;x_3) = 0,693147 + (x - 2) \cdot 0,405466+\\ + (x-2)(x-3) \cdot (-0,058892) + (x-2)(x-3)(x-4) \cdot 0,0887416.

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

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

 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

Страницы: 1, 2


© 2010 САЙТ РЕФЕРАТОВ