">
Математика Математический анализ
Информация о работе

Тема: Решение краевой задачи дляобыкновенного дифференциального уравнения.

Описание: Постановка задачи. Метод сведения краевой задачи к задаче Коши, конечных разностей. Описание его, результатов. Метод Галёркина, описание. Листинг программы. Линейное дифференциальное уравнение. Проверка достигнутой точности. Соответствующие значения функции.
Предмет: Математика.
Дисциплина: Математический анализ.
Тип: Курсовая работа
Дата: 31.08.2012 г.
Язык: Русский
Скачиваний: 3
Поднять уникальность

Похожие работы:

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

Решение краевой задачи для

обыкновенного дифференциального уравнения.

Вариант №1

2008

Содержание:

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

Метод сведения краевой задачи к задаче Коши:

2.1. Описание метода..............................................................................................................4

2.2. Описание результатов………………………………………….....................................7

Метод конечных разностей:

3.1. Описание метода …………………………………………………………….............…8

3.2. Описание результатов...................................................................................................11

Метод Галёркина:

4.1. Описание метода ………………………………….…………………………………..12

4.2. Описание результатов...................................................................................................15

Выводы..………………..…………………………………………………………………...16

Приложение 1. Листинг программы «Метод сведения краевой задачи к задаче Коши»……..17

Приложение 2. Листинг программы «Метод конечных разностей»…………………………...19

Приложение 3. Листинг программы «Метод Галёркина»………………………………………21

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

Дано линейное дифференциальное уравнение второго порядка:

 (1.1)

с краевыми условиями:

(1.2)

где функции P(x), Q(x), F(x) непрерывны, и - заданные постоянные, причем:  и .

Требуется найти решение y(x) уравнения (1.1), удовлетворяющее краевым условиям (1.2).

В данном варианте курсовой работы необходимо:

Решить краевую задачу для обыкновенного дифференциального уравнения с точностью Е=1.е-4. Проверить достигнутую точность. Результаты представить с шагом h=0.04. Решение провести тремя следующими методами:

1. Сведением краевой задачи к задаче Коши.

Методом конечных разностей.

Методом Галёркина.

 (1.3)

Из условий следует, что функции P(x) = -4/x; Q(x) =-4/x ; F(x) =-1/x-2.8 и постоянные

(1 = 1, (2 = 0, (1 = 1, (2 = 0, a =2.6, b =3.4, A = -3.8, B = 1.

2. Метод сведения краевой задачи к задаче Коши

2.1. Описание метода

Решение дифференциального уравнения (1.1) с краевыми условиями (1.2) будем искать в виде линейной комбинации:

где c=const(2.1)

Подставим y(x) в виде (2.1) в исходное дифференциальное уравнение (1.1) и сгруппируем слагаемые при постоянной с, получим выражение:

(2.2)

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

(2.3)

Из системы дифференциальных уравнений (2.3) видно, что функция u=u(x) - ненулевое решение соответствующего однородного уравнения, а v=v(x) - некоторое решение данного неоднородного уравнения (1.1).

Чтобы свести краевую задачу (1.1)-(1.2) к задачам Коши для функций u=u(x) и v=v(x), подставим в первое краевое условие (1.2) выражение для функции y(x) и сгруппируем слагаемые при постоянной c, будем иметь:

 (2.4)

Для того чтобы равенство (2.4) было справедливо при любом с, необходимо и достаточно, чтобы коэффициенты при постоянной с обращались в ноль, т. е. должны быть выполнены равенства:

 (2.5)

Получили систему (2.5) из двух уравнений с четырьмя неизвестными. Решений системы будет бесконечное множество. Найдем хотя бы одно.

Для обеспечения первого равенства системы, например, можно подобрать:

 (2.6)

где постоянная k отлична от нуля, так как тривиальное решение u(a)=0 можно отбросить.

Для выполнения второго равенства системы (2.5) можно положить

   (2.7)

или

при (2.8)

По условию курсовой работы (1 = 1, (2 = 0, A = -3.8 будем иметь:

(2.9)

 (2.10)

Отсюда видно, что u есть решение задачи Коши (2.9) для однородного уравнения (2.3), удовлетворяющее начальным условиям (2.6), а v есть решение задачи Коши (2.10) для неоднородного уравнения (2.3), удовлетворяющее начальным условиям (2.7) или (2.8). При этом для любого с функция y = cu + v удовлетворяет краевому условию на конце x=a.

Подберем теперь постоянную c так, чтобы функция y(x) удовлетворяла краевому условию (1.2) на конце x=b. Это дает:



откуда:

,

при этом предполагается, что знаменатель

Так как по условию работы (1 = 1, (2 = 0, то коэффициент c будет иметь вид:



Для решения полученных уравнений из системы (2.3) будем использовать метод Рунге-Кутта, который имеет достаточно высокую точность на всем интервале порядка 

Рассмотрим второе дифференциальное уравнение из системы (2.3) с начальным условием (2.7). Дифференциальные уравнения системы (2.3) являются однотипными, поэтому решение первого уравнения системы (2.3) с начальным условием (2.6) осуществляется аналогичным способом, при условии .

Для того чтобы решить указанное дифференциальное уравнение методом Рунге-Кутта, сделаем замену:

 

и подставим ее во второе дифференциальное уравнение из системы (2.3), получим:

, 

Расчет производим следующим образом: выберем шаг h, приращение x в зависимости от шага будет:



где n=0,1…k-1

Соответствующие значения  и искомых функций v и z определяются формулами:



где:





Чтобы достичь заданную точность вычисляем y(x) двумя способами: один раз с шагом h, другой раз с шагом h/2, получая при этом значения более точные. Если расхождение полученных значений не превышает заданной точности Е =, то выбранный шаг h можно считать достаточным и полученная функция y(x) удовлетворяет заданной точности. Иначе уменьшаем шаг h, пока не будет достигнута заданная точность.

2.1. Описание результатов

При решении данного дифференциального уравнения второго порядка с заданными краевыми условиями (1.3) методом сведения к задаче Коши и последующим её решением методом Рунге-Кутта, были получены следующие результаты, представленные в таблице 1. В столбце Х приведено разбиение отрезка [2.6; 3.4] с шагом h = 0.04, в столбце Y(X) – значение функции  (n=0,…,20) в соответствующих точках , в столбце Delta – значения найденных абсолютных погрешностей.

В результате работы программы, листинг которой приведен в приложении 1, точность была достигнута при шаге 0.02, количество разбиений отрезка [2.6; 3.4] равно 40, максимальная погрешность равна как видно из таблицы 1 при n=11 или 12.

Точность решений определялась из условия, что норма модуля разности более точного решения и приближенного должна быть меньше заданной точности. Для достижения заданной точности E =  шаг h = 0.04 уменьшили в 2 раза, так как при шаге h = 0.04 полученная точность не удовлетворяла заданной.

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

Таблица 1

n X Y(X) Delta  0 2.60 -3.800000 0  1 2.64 -3.582867   2 2.68 -3.366368   3 2.72 -3.150070   4 2.76 -2.933525   5 2.80 -2.716271   6 2.84 -2.497832   7 2.88 -2.277714   8 2.92 -2.055408   9 2.96 -1.830382   10 3.00 -1.602085   11 3.04 -1.369944   12 3.08 -1.133363   13 3.12 -0.891717   14 3.16 -0.644358   15 3.20 -0.390608   16 3.24 -0.129757   17 3.28 0.138936   18 3.32 0.416246   19 3.36 0.702984   20 3.40 1.000000 0  

3. Метод конечных разностей

3.1. Описание метода

Решим дифференциальное уравнение (1.1) с краевыми условиями (1.2) методом конечных разностей. Для того чтобы получить систему конечно-разностных уравнений, разобьем отрезок [a,b] на n равных частей длины  как показано на рисунке 1. Точки разбиения имеют абсциссы:

 (i=0,1,2,…,n),  .



Рис. 1

Значения в точках  искомой функции y=y(x) и ее производных обозначим соответственно:



Введем также обозначения: .

Заменяя производные симметричными конечно-разностными отношениями для внутренних точек  отрезка [a, b] будем иметь:

 (i=1,2,…,n-1) (3.1)

где - погрешность формулы порядка .

Подставим конечно-разностные выражения (3.1) в исходное дифференциальное уравнение (1.1). В результате подстановки получим:

 (3.2)

Или запишем выражение (3.2) в виде:

,(i=1,2,...,n-1), (3.3)

где  

Так как всего n+1 неизвестных  на отрезке [a, b], а уравнений n-1, то не хватает еще двух уравнений для нахождения всех  (i=0,1,…,n). Эти недостающие уравнения дают краевые условия (1.2). По условию курсовой работы , , поэтому система (1.2) принимает вид:



С учетом этого запишем систему линейных алгебраических уравнений:

, где (i=1,2,...,n-1) (3.4)

Таким образом, получена линейная система (3.4) из n+1 уравнений с n+1 неизвестными , представляющими собой значения искомой функции y=y(x) в точках .

Систему уравнений (3.4) можно записать в матричном виде, при этом получается трехдиагональная матрица коэффициентов:

1 0 0 0 0 ..…………. 0 0 0 0  A

  1 0 0 ………….....0 0 0 0 

0   1 0………….....0 0 0 0 ……. …….

………………………….……....……... …….. = …….. (3.5)

………………………….……….....….. …….. …….

………………………….......1 ……. 

…………………………….…...... 0 1  B

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

Приведем матрицу (3.5) к двухдиагональному виду. Предположим, что из (3.4) исключена неизвестная . Тогда это уравнение примет вид:

 (3.7)

где  (i=1,2,...,n-1) некоторые коэффициенты.

В прямом ходе метода прогонки определяются коэффициенты  зная  Найдем эти коэффициенты, для этого запишем



Подставим это выражение в (3.4) и, выразив , будем иметь:

 (3.8)

Сравнивая формулы (3.7) и (3.8), получим рекуррентные формулы для определения  и :

 (i=1,2,...,n-1)(3.9)

Определим теперь  и . Из (3.7) при i=0 получим:

 (3.10)

Сравнивая (3.10) и первое краевое условие из (1.2), находим:



При обратном ходе, используя формулу (3.7) и условие , последовательно находим 

Точность в методе конечных разностей достигается аналогичным способом, приведенным в методе сведения краевой задачи к задаче Коши, так как полученная точность зависит от шага разбиения отрезка [a, b].

3.2. Описание результатов

При решении данного дифференциального уравнения второго порядка с заданными краевыми условиями (1.3) методом конечных разностей, получены следующие результаты представленные в таблице 2. В столбце Х приведено разбиение отрезка [2.6; 3.4] с шагом h = 0.04, в столбце Y(X) – значение функции  (n=0,…,20) в соответствующих точках , в столбце Delta - значения найденных абсолютных погрешностей.

В результате работы программы, листинг которой приведен в приложении 2, точность была достигнута при шаге 0.01, количество разбиений отрезка [2.6; 3.4] равно 80, максимальная погрешность как видно из таблицы 2 при n=10,11 или 12.

Для достижения заданной точности Е = шаг h = 0.04 уменьшили в 4 раза, поскольку при шаге большем, чем h = 0.01 полученная точность не удовлетворяла заданной.

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

Таблица 2 n X Y(X) Delta  0 2.60 -3.800000 0  1 2.64 -3.582868   2 2.68 -3.366371   3 2.72 -3.150074   4 2.76 -2.933530   5 2.80 -2.716277   6 2.84 -2.497839   7 2.88 -2.277722   8 2.92 -2.055416   9 2.96 -1.830391   10 3.00 -1.602094   11 3.04 -1.369954   12 3.08 -1.133372   13 3.12 -0.891726   14 3.16 -0.644366   15 3.20 -0.390615    16 3.24 -0.129763   17 3.28 0.138931   18 3.32 0.416242   19 3.36 0.702982   20 3.40 1.000000 0  

4. Метод Галёркина

4.1. Описание метода

В данном методе решения дифференциального уравнения (1.1) с краевыми условиями (1.2) для удобства описания метода введем:

линейный дифференциальный оператор ,

линейные операторы 

Тогда постановка задачи для дифференциального уравнения (1.1) с краевыми условиями (1.2) примет вид:

 (4.1)

Решение краевой задачи (4.1) будем искать в виде суммы:

 (4.2)

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

Функция y(x) должна удовлетворять краевым условиям, то есть должно выполняться:



Краевые условия должны быть справедливы для любого набора констант , поэтому запишем:

 

По условию курсовой работы  получим:

  (4.3)

В качестве функций  будем выбирать функции, принадлежащие к классу полиномов. Тогда, учитывая краевые условия, запишем полином Лагранжа для  получим:



В качестве базисных функций  выбираем:



При таком подборе базисных функций  функция y(x), удовлетворяет краевым условиям, при любом наборе констант .

Введем следующее определение: невязка – разность между левой и правой частями уравнения:



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

Таким образом, для определения коэффициентов  (i=1,2,...,n) приходим к системе линейных уравнений:

, (i=1,2,...,n)

или более подробно запишем данную систему в матричном виде:

 = (4.4)

Система уравнений (4.4) разрешается относительно  по методу Гаусса. По полученным значениям  и функциям  определяем столбец значений функции y(x) согласно формуле (4.2).

Определенные интегралы системы уравнений (4.4) можно вычислить по методу Симпсона, общая формула которого имеет вид:



где n - количество разбиений отрезка [a, b].

Общая погрешность метода Симпсона составляет  где .

Подынтегральные функции в системе (4.4) имеют видили  Тогда из однородных краевых условий (4.3) получим  или , что упрощает формулу Симпсона.

Чтобы достичь заданную точность вычисляем y(x) двумя способами: один раз с числом базисных функций i, другой раз с i+1, получая при этом значения более точные. Если расхождение полученных значений не превышает заданной точности Е = то выбранное число функций можно считать достаточным и полученная функция y(x) удовлетворяет заданной точности. Иначе увеличиваем i, пока не будет достигнута заданная точность.

4.1. Описание результатов

При решении данного дифференциального уравнения второго порядка с заданными краевыми условиями (1.3) методом Галёркина, получены следующие результаты представленные в таблице 3. В столбце Х приведено разбиение отрезка [2.6; 3.4] с шагом h = 0.04, в столбце Y(X) – значение функции  (n=0,…,20) в соответствующих точках в столбце Delta – значения найденных абсолютных погрешностей.

В результате работы программы, листинг которой приведен в приложении 3, точность была достигнута при количестве базисных функций равным 5, максимальная погрешность  как видно из таблицы 3 при n=1.

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

Таблица 3

n X Y(X) Delta  0 2.60 -3.800000 0  1 2.64 -3.582867   2 2.68 -3.366369   3 2.72 -3.150070   4 2.76 -2.933525   5 2.80 -2.716271   6 2.84 -2.497832   7 2.88 -2.277714   8 2.92 -2.055408   9 2.96 -1.830382   10 3.00 -1.602085   11 3.04 -1.369944   12 3.08 -1.133363    13 3.12 -0.891717   14 3.16 -0.644358   15 3.20 -0.390608   16 3.24 -0.129757   17 3.28 0.138936   18 3.32 0.416246   19 3.36 0.702984   20 3.40 1.000000 0  

5. Выводы

Обыкновенное дифференциальное уравнение второго порядка было решено тремя методами: сведением краевой задачи к задаче Коши, методом конечных разностей, методом Галёркина. Для получения решения были написаны программы на языке программирования C++.

Максимальная погрешность в методе сведения краевой задачи к задаче Коши равна , в методе конечных разностей: , в методе Галёркина: .

Требуемая точность E=была достигнута. При использовании методов сведения к задаче Коши и конечных разностей точность достигалась путём уменьшения шага, при использовании метода Галёркина - путём увеличения количества базисных функций.

Приложение 1.

Листинг программы «Метод сведения краевой задачи к задаче Коши»

#include "stdio.h"

#include "math.h"

const int N_max=100;

const double a=2.6, b=3.4, A=-3.8, B=1, epsel=0.0001;

double F(double x)

{ return -1/x-2.8;}

double Q(double x)

{ return -4/x;}

double P(double x)

{ return -4/x;}

void Runge_Kutt(double h,int n,double U[],double V[],double z[],double t[],double a)

{

double k1z,k2z,k3z,k4z,

k1v,k2v,k3v,k4v,

k1t,k2t,k3t,k4t,

k1u,k2u,k3u,k4u;

double x=a;

z[0]=0;

V[0]=A;

t[0]=1;

U[0]=0;

for(int i=0;i<=n;i++){

x=a+i*h;

k1v=z[i];

k1z=F(x) - P(x)*z[i]-Q(x)*V[i];

k2v=z[i]+h/2*k1z;

k2z=F(x+h/2) - P(x+h/2)*(z[i]+h/2*k1z)-Q(x+h/2)*(V[i]+h/2*k1v);

k3v=z[i]+h/2*k2z;

k3z=F(x+h/2) - P(x+h/2)*(z[i]+h/2*k2z)-Q(x+h/2)*(V[i]+h/2*k2v);

k4v=z[i]+h*k3z;

k4z=F(x+h) - P(x+h)*(z[i]+h*k3z)-Q(x+h)*(V[i]+h*k3v);

k1u=t[i];

k1t=-P(x)*t[i]-Q(x)*U[i];

k2u=t[i]+h/2*k1t;

k2t=-P(x+h/2)*(t[i]+h/2*k1t)-Q(x+h/2)*(U[i]+h/2*k1u);

k3u=t[i]+h/2*k2t;

k3t=-P(x+h/2)*(t[i]+h/2*k2t)-Q(x+h/2)*(U[i]+h/2*k2u);

k4u=t[i]+h*k3t;

k4t=-P(x+h)*(t[i]+h*k3t)-Q(x+h)*(U[i]+h*k3u);

V[i+1]=V[i]+h/6*(k1v+2*k2v+2*k3v+k4v);

z[i+1]=z[i]+h/6*(k1z+2*k2z+2*k3z+k4z);

U[i+1]=U[i]+h/6*(k1u+2*k2u+2*k3u+k4u);

t[i+1]=t[i]+h/6*(k1t+2*k2t+2*k3t+k4t);

}

}

double Proverka_Tochnosti(double y[],double y_tochnoe[],int n,double delta,double eps[])

{

int i=0;

delta=0;

for(i=0;i<=n/2;i++)

{

eps[i]=fabs(y[2*i]-y_tochnoe[i]);

if(eps[i]>delta)

delta=eps[i];

}

return delta;

}

void Nahozdenie_y(double y[],double U[],double V[],int n,double C)

{

int i=0;

for(i=0;i<=n;i++)

{y[i]=C*U[i]+V[i];}

}

void Prisvoenie_y_Tochnoe(double y[],double y_tochnoe[],int n)

{

int i=0;

for(i=0;i<=n;i++)

{y_tochnoe[i] = y[i];}

}

void Pechat(double y_tochnoe[],double a,double eps[],int n,int N,double h,double h1,double delta)

{

int i=0,j=0;

double x=a;

printf("n = %d h = %1.2f Max_Delta = %1.2e ",n/2,h1,delta);

printf("XY(X)Delta ");

for(i=0;i<=n/2;i=i+n/2/N)

{

x=a+j*h;

j++;

printf("%1.2f%f%1.2e ",x,y_tochnoe[i],eps[i/2]);

}

getchar();

}

main()

{

int i=0,j=0,N=20,n=0;

double U[N_max],V[N_max],z[N_max],y[N_max],y_tochnoe[N_max],t[N_max],eps[N_max],C=0;

double h=0,x=0,delta=0;

n=N;

do

{

if(n>=N_max)

{

printf("Error!!! Vihod za predeli massiva");

return 0;

}

h=(b-a)/n;

Runge_Kutt(h,n,U,V,z,t,a);

C=(B-V[n])/U[n];

Nahozdenie_y(y,U,V,n,C);

delta=Proverka_Tochnosti(y,y_tochnoe,n,delta,eps);

Prisvoenie_y_Tochnoe(y,y_tochnoe,n);

n=2*n;

}while(delta>epsel);

Pechat(y_tochnoe,a,eps,n,N,0.04,h,delta);

return 0;}

Приложение 2.

Листинг программы «Метод конечных разностей»

#include "stdio.h"

#include "math.h"

const int N_max=200;

const double epsel=0.0001;

const double a=2.6, b=3.4, A=-3.8, B=1;

double F(double x)

{ return -1/(x)-2.8;}

double Q(double x)

{ return -4/(x);}

double P(double x)

{ return -4/(x);}

void Prisvoenie_y_Tochnoe(double y[],double y_tochnoe[],int n)

{

int i=0;

for(i=0;i<=n;i++)

{y_tochnoe[i] = y[i];}

}

void Konec_Raznost(double h,double m[],double k[],double f[],double c[],double d[],double temp[],int n )

{

double x=0;

c[0]=0;

temp[0]=A;

for(int i=0;i {

x=a+i*h;

m[i]=(Q(x) - 2/(h*h)) / (1/(h*h) + P(x)/(2*h));

k[i]= (1/(h*h) - P(x)/(2*h)) / (1/(h*h) + P(x)/(2*h));

f[i]=F(x) / (1/(h*h) + P(x)/(2*h));

}

for(int i=1;i {

c[i]=1/(m[i]-k[i]*c[i-1]);

d[i]=f[i]-k[i]*temp[i-1];

temp[i]=c[i]*d[i];

}

}

double Proverka_Tochnosti(double y[],double y_tochnoe[],int n,double delta,double eps[])

{

int i=0;

delta=0;

for(i=0;i<=n/2;i++)

{

eps[i]=fabs(y[2*i]-y_tochnoe[i]);

if(eps[i]>delta)

delta=eps[i];

}

return delta;

}

void Nahozdenie_y(double y[],double c[],double temp[],int n)

{

int i=0;

y[n]=B;

for(i=n;i>0;i--)

{y[i-1]=temp[i-1]-c[i-1]*y[i];}

}

void Pechat(double y_tochnoe[],double a,double eps[],int n,int N,double h,double h1,double delta)

{

int i=0,j=0;

double x=a;

printf("n = %d h = %1.2f Max_Delta = %1.2e ",n/2,h1,delta);

printf("XY(X)Delta ");

for(i=0;i<=n/2;i=i+n/2/N)

{

x=a+j*h;

j++;

printf("%1.2f%f%1.2e ",x,y_tochnoe[i],eps[i/2]);

}

getchar();

}

main()

{

int i=0,j=0,N=20,n=0;

double y[N_max],y_tochnoe[N_max],eps[N_max],C=0;

double h=0,x=0,delta=0;

double m[N_max],k[N_max], f[N_max], c[N_max], d[N_max], temp[N_max];

n=N;

do

{

if(n>=N_max)

{

printf("Error!!! Vihod za predeli");

return 0;

}

h=(b-a)/n;

Konec_Raznost(h,m,k,f,c,d,temp,n);

Nahozdenie_y(y,c,temp,n);

delta=Proverka_Tochnosti(y,y_tochnoe,n,delta,eps);

Prisvoenie_y_Tochnoe(y,y_tochnoe,n);

n=2*n;

}while(delta>epsel);

Pechat(y_tochnoe,a,eps,n,N,0.04,h,delta);

return 0;}

Приложение 3.

Листинг программы «Метод Галёркина»

#include "stdio.h"

#include "math.h"

#include "process.h"

#include "conio.h"

const int N_max=21,E=8;

const double epsel=0.0001,h=0.04;

const double a=2.6, b=3.4, A=-3.8, B=1;

void Prisvoenie_y_tochnoe(double y[],double y_tochnoe[],int n)

{

int i=0;

for(i=0;i {y_tochnoe[i] = y[i];}

}

double F(double x)

{return -1/(x)-2.8;}

double Q(double x)

{return -4/(x);}

double P(double x)

{return -4/(x);}

double Proverka(double y[],double y_tochnoe[],int n,double delta,double eps[])

{

int i=0;

delta=0;

for(i=0;i {

eps[i]=fabs(y[i]-y_tochnoe[i]);

if(eps[i]>delta)

delta=eps[i];

}

return delta;

}

void Nahozdenie_y(double y[],double C[],int n)

{

int i=0;

double x=0;

for(int j=0;j x=a+j*h;

y[j]=(A*b-B*a)/(b-a)+(B-A)/(b-a)*x;

for(i=1;i<=n;i++)

{y[j]+=C[i-1]*pow(x-a,i)*(x-b);}

}

}

void Pechat_rez(double y_tochnoe[],double eps[], int n,double delta)

{

int i=0,j=0;

double x=a;

printf("Bazis = %d ",n);

printf("Razbienij v metode Simpsona = %d ",E);

printf("Max_Delta = %1.2e ",delta);

printf(" X Y(X) Delta ");

for(j=0;j {

x=a+j*h;

printf(" %1.2f %lf %1.2e ",x,y_tochnoe[j],eps[j]);

}

getchar();

}

void Resheniya(double matrix[N_max][N_max],double D[],double C[],int n)

{

double U,L;

double x=0,sum1=0,sum2=0,h=0;

int k=0,i=0,j=0;

h=(b-a)/E;

for(k=1;k<=n;k++){

for(i=1;i<=n;i++)

{

for(j=1;j {

x=a+j*h;

U=pow(x-a,i)*(x-b);

L=k*pow(x-a,k-2)*((k-1)*(x-b)+2*(x-a))+P(x)*(k*pow(x-a,k-1)*(x-b)+pow(x-a,k))+Q(x)*pow(x-a,k)*(x-b);

if(j%2) sum1+=U*L;

else sum2+=U*L;

}

matrix[i-1][k-1]=h/3*(4*sum1+2*sum2);

sum1=0;sum2=0;

}

}

for(i=1;i<=n;i++)

{

for(j=1;j {

x=a+j*h;

U=pow(x-a,i)*(x-b);

L=P(x)*(B-A)/(b-a)+Q(x)*((A*b-B*a)/(b-a)+(B-A)/(b-a)*x);

if(j%2) sum1+=U*(F(x)-L);

else sum2+=U*(F(x)-L);

}

D[i-1]=h/3*(4*sum1+2*sum2);

sum1=0;sum2=0;

}

// Gauss

double koef=0;

for(int t=0;t {

for( i=t+1;i koef=0;

koef=matrix[i][t]/matrix[t][t];

D[i]=D[i]-D[t]*koef;

for( j=t+1;j matrix[i][j]=matrix[i][j]-matrix[t][j]*koef;

}

}

}

for(i=n-1;i>=0;i--){

C[i]=D[i];

for(j=i+1;j C[i]=C[i]- matrix[i][j]*C[j];

}

C[i]=C[i]/matrix[i][i];

}

}

main()

{

int i=0,j=0,n=0;

double y[N_max],y_tochnoe[N_max],eps[N_max];

double matrix[N_max][N_max],B[N_max],C[N_max];

double x=0,delta=0;

n=0;

do

{

n++;

if(n>=10)

{

printf("Error!!! Vihod za predeli ");

getchar();

return 0;

}

double x=0,sum1=0,sum2=0;

int k=0,i=0,j=0;

Resheniya(matrix,B,C,n);

Nahozdenie_y(y,C,n);

delta=Proverka(y,y_tochnoe,n,delta,eps);

Prisvoenie_y_tochnoe(y,y_tochnoe,n);

}while(delta>epsel);

Pechat_rez(y_tochnoe,eps,n,delta);

return 0;

}