Решение краевой задачи методом прогонки (ПЗ + программа на С++)

- Москва : Наука, 1978.


Приложение
//Описание используемых библиотек
//Описание билиотек ввода вывода
#include
#include
//Опиасание библиотеки математических функций
#include
//Описание стандартных библиотек
#include
#include
//Описание области видимости функций из библиотеки stdlib
using namespace std;
//функция считывания параметров задачи
void read(double &a, double &b, int &n, double &alfa0, double &alfa1, double &A1, double &betta0, double &betta1, double &B1)
{
//открываем входной поток
ifstream f("data.txt");
//считываем левую границу
f >> a;
//считываем правую границу
f >> b;
//считываем число интервалов разбиения
f >> n;
//считываем параметры левого краевого условия
f >> alfa0;
f >> alfa1;
f >> A1;
//считываем параметры правого краевого условия
f >> betta0;
f >> betta1;
f >> B1;
//закрываем входной поток
f.close();
};
//Функция считывания значений табличной функции p(x)
void read_p(int n, int j, double **f) //n-количество точек,
//f-массив значений функции
//j-количество вариантов функции (число столбцов)
{
ifstream f1("f.txt"); //открываем входной поток
int i, k; //счетчики
//считываем табличную функцию
for(i = 0; i < n; i++)
for(k = 0; k < j; k++)
f1 >> f[i][k];
//Закрываем входной поток
f1.close();
};
//функция вычисления q(x)
void q_rechnen(double a, double b, int n, double u, double *q)
//a-левая граница, b-правая граница, n-число интервалов,
//u-параметр функции, q-массив значений функции в точках
{
int i; //счетчик
double h = (b - a)/n; //шаг интрвалов
for(i = 0; i <= n; i++)
q[i] = cos(u * (a + i * h));
};
//функция вычисления правой части f(x)
void f_rechnen(double a, double b, int n, double *f)
//a-левая граница, b-правая граница, n-число интервалов,
//f-массив значений функции в точках
{
int i; //счетчик
double h = (b - a)/n; //шаг интрвалов
double x; //координаты узлов
for(i = 0; i <= n; i++)
{
x = a + h * i;
f[i] = 2 * x * x + 2 * x - 4;
}
};
//головная функция
void main()
{
double a, b; //левая правая граница
int n; //число интервалов
double alfa0, alfa1, A1; //параметры левого краевого условия
double betta0, betta1, B1; //параметры правого краевого условия
//чтение параметров задачи
read(a, b, n, alfa0, alfa1, A1, betta0, betta1, B1);
double h = (b - a)/n; //вычисление шага
int j = 5, s, k, i; //описание счетчиков
double u; //параметр задачи
//описание массивов
double *x; //массив узлов
x = new double[n + 1];
//вычисление массива узлов
for(i = 0; i <= n; i++)
x[i] = a + i * h;
//массив значений функции p(x) табличная
double **p;
p = new double*[n + 1];
for(i = 0; i < n + 1; i++)
p[i] = new double[j];
read_p(n, j, p);
//массив значений правой части
double *f;
f = new double[n + 1];
f_rechnen(a, b, n, f);
//массив значений функции q(x)
double *q;
q = new double[n + 1];
//описание коэффициентов прогонки
double *An, *Bn;
An = new double[n];
Bn = new double[n];
//описание массива решений
double *y;
y = new double[n + 1];
//описание массивов матрицы
double *c_n, *a_n, *b_n, *d_n;
c_n = new double[n + 1];
a_n = new double[n + 1];
b_n = new double[n + 1];
d_n = new double[n + 1];
//отрытие выходного потока
ofstream f1("result.txt");
for(s = 0; s < 6; s++) //счетчик по параметру s
{
u = 0.7 + 0.05 * s; //вычислим параметр u
//вычислим значение функции q(u,x) в узлах сетки
q_rechnen(a, b, n, u, q);
//Вычислим элементы СЛАУ на краевых условиях:
b_n[0] = h * alfa0 - alfa1;
c_n[0] = alfa1;
d_n[0] = A1 * h;
a_n[0] = 0.0;
a_n[n] = -betta1;
b_n[n] = h * betta0 + betta1;
c_n[n] = 0.0;
d_n[n] = h * B1;
for(k = 0; k < j; k++) //цикл по функциям p_j
{
//вычисляем внекраевые элементы матрицы СЛАУ
for(i = 0; i < n - 2; i++)
{
a_n[i + 1] = 1.0 - h * p[i][k] + h * h * q[i];
b_n[i + 1] = -2.0 + h * p[i][k];
c_n[i + 1] = 1;
d_n[i + 1] = h * h * f[i];
}
//вычисляем коэффициенты прогонки
An[0] = -alfa1/(h * alfa0 - alfa1);
Bn[0] = A1 * h/(h * alfa0 - alfa1);
for(i = 1; i < n; i++)
{
An[i] = -1.0/(a_n[i] * An[i - 1] + b_n[i]);
Bn[i] = (d_n[i] - a_n[i] * Bn[i - 1])/(a_n[i] * An[i - 1] + b_n[i]);
}
//Вычисляем значение функции в последней точке;
y[n] = (h * B1 + betta1 * Bn[n - 1])/(h * betta0 + betta1 * (1 - An[n - 1]));
for(i = n - 1;i >= 0; i--)
y[i] = An[i] * y[i + 1] + Bn[i];
//печатаем результат решения в файл
f1 <<"j="<< k + 1 <<" s="<< s << endl;
for(i = 0; i <= n; i++)
f1 << x[i] <<"\t" << y[i] << endl;
}
}
};








7