Задача 5.2
Постановка задачи
Функцияu(x), являющаяся решением краевой задачи
описывает стационарное распределение тепла в одномерном стержне, совпадающем с отрезком 
.
Рассчитать стационарное поле температур на отрезке . Предусмотреть возможность задания произвольных граничных условий. Для проверки работы программы построить и рассчитать тестовый пример.
| 5.2.11 | 2.2 | 
| 
|
Решение
1. Построим тестовый пример.
Выберем в качестве точного решения некоторую непрерывную функцию: u(x)=1-x2. Очевидно, что значение u(a)=u(0)= -1. Это значит, что в тестовом примере ua=1. Аналогично, ub=-3.84.
Построим правую часть тестового примера. Для этого найдем левую часть уравнения при выбранной функции u(x).
g(x) = -u”(x) + q(x)*u(x) = 2 + (2+x)^(1/2)*(1-x2) = 2+(2+x)^(1/2)-(x^2)* *(2+x)^(1/2).
Этому же должна быть равна и правая часть уравнения, составляющая сейчас:
f(x) = (u(x)^2 + u(x)) = 2 – 3*x^2 + x^4.
Ясно, что равенства можно добиться, добавляя в правую часть слагаемое
s(x) = g(x) – f(x) =(2+x)^(1/2) - (x^2)*(2+x)^(1/2) + 3*x^2 - x^4
Таким образом, построен тестовый пример:
-u”+xu = u^2+u + (2+x)^(1/2) - (x^2)*(2+x)^(1/2) + 3*x^2 - x^4 x ϵ (0, 2.2),
u(0) = 1,
u(2.3) = -3.84.
Решением построенной задачи является функция u(x)=1-x2.
2. Проверим правильность работы программы на тестовом примере.
Посчитаем значения функции u(x) в точках 0.575, 1.15, 1.725.
u(0.55) = 1-0.552 = 0.6975,
u(1.1) = 1 - 1.12 = -0.21,
u(1.65) = 1 - 1.652 = -1.7225.
Результаты программы совпадают с вычислениями, значит, программа работает корректно.
Введите количество разбиений отрезка [a,b]
Введите точность эпсилон
0.00001
Введите значение в точке а
Введите значение в точке b
-3.84
Ua: 1, Ub: -3.84
Деления отрезка:
0 0.55 1.1 1.65 2.2
Матрица производных:
1 0 0 0 0
-1 2.78555 -1 0 0
0 -1 2.83511 -1 0
0 0 -1 2.88043 -1
0 0 0 0 1
Наддиагональ:
0 -1 -1 -1 0
Поддиагональ:
0 -1 -1 -1 0
Главная диагональ:
1 2.78555 2.83511 2.88043 1
Правая часть:
0 3.06682 1.07594 -3.02901 0
Начальное значение:
1 -1 -1 -1 -3.84
Ответ:
1 0.6975 -0.21 -1.7225 -3.84
Для продолжения нажмите любую клавишу...
3. Решим данную задачу. Граничные условия зададим такие: u(0) = -1, u(2.2) = -11.
Введите количество разбиений отрезка [a,b]
Введите точность эпсилон
0.000001
Введите значение в точке а
-1
Введите значение в точке b
-11
Ua: -1, Ub: -11
Деления отрезка:
0 0.244444 0.488889 0.733333 0.977778 1.22222 1.46667 1.71111 1.95556 2.2
Матрица производных:
1 0 0 0 0 0 0 0 0 0
-1 2.14927 -1 0 0 0 0 0 0 0
0 -1 2.15402 -1 0 0 0 0 0 0
0 0 -1 2.15854 -1 0 0 0 0 0
0 0 0 -1 2.16286 -1 0 0 0 0
0 0 0 0 -1 2.16701 -1 0 0 0
0 0 0 0 0 -1 2.17101 -1 0 0
0 0 0 0 0 0 -1 2.17486 -1 0
0 0 0 0 0 0 0 -1 2.17859 -1
0 0 0 0 0 0 0 0 0 1
Наддиагональ:
0 -1 -1 -1 -1 -1 -1 -1 -1 0
Поддиагональ:
0 -1 -1 -1 -1 -1 -1 -1 -1 0
Главная диагональ:
1 2.14927 2.15402 2.15854 2.16286 2.16701 2.17101 2.17486 2.17859 1
Правая часть:
0 0.0895189 0.0942677 0.0987886 0.103111 0.10726 0.111254 0.11511 -9.88116 0
Начальное значение:
-1 -1 -1 -1 -1 -1 -1 -1 -1 -11
Ответ:
-1 -0.958525 -1.00048 -1.13678 -1.39467 -1.82925 -2.55068 -3.79223 -6.10301 -11
Для продолжения нажмите любую клавишу...
4. Листинг программы (только существенные части):
// Значения, стоящие в матрице производных на главной диагонали
void diagElem(double mas[], double z[], double x[], int N)
{
double a = 0; // левая граница стержня
double b = 2.2; // правая граница стержня
double h = (b - a) / N; // шаг
for (int i = 0; i <= N; i++)
{
mas[0] = 1;
if ((i!= 0) & (i!= N))
mas[i] = (2 + h*h*(sqrt(2 + x[i]))) - h*h*(2 * z[i] + 1);
mas[N] = 1;
}
}
// Функция (правая часть м.Ньютона)
void F(double mas[], double z[], double x[], int N, double grA, double grB)
{
double a = 0; // левая граница стержня
double b = 2.2; // правая граница стержня
double h = (b - a) / N; // шаг
mas[0] = -z[0] + grA;
for (int i = 1; i < N; i++)
{
mas[i] = -(-z[i - 1] + 2 * z[i] - z[i + 1] + sqrt(2 + x[i])*z[i] * h*h - h*h*(z[i] * z[i] + z[i]));
}
mas[N] = -z[N] + grB;
}
// Метод прогонки
void solveProgonka(double x[], double a[], double b[], double c[], double d[], int n)
{
double alpha[nmax], betta[nmax];
alpha[0] = -c[0] / b[0];
betta[0] = d[0] / b[0];
for (int i = 1; i < n; i++)
{
alpha[i] = -c[i] / (b[i] + a[i] * alpha[i - 1]);
betta[i] = (d[i] - a[i] * betta[i - 1]) / (b[i] + a[i] * alpha[i - 1]);
}
betta[n] = (d[n] - a[n] * betta[n - 1]) / b[n] + a[n] * alpha[n - 1];
x[n] = betta[n];
for (int i = n - 1; i >= 0; i--)
{
x[i] = alpha[i] * x[i + 1] + betta[i];
}
}
// Первая норма вектора
double norm1(double x[], int n)
{
double sum;
sum = 0.0;
for (int i = 0; i <= n; i++)
sum = sum + abs(x[i]);
return sum;
}
// Цикл поиска решения
while (norma>eps)
{
// сохраняем старое значение
for (inti = 0; i<= N; i++)
temp[i] = z[i];
// изменяем матрицу производных
diagElem(mainDiag, z, setka, N);
// изменяем правую часть
F(funcF, z, setka, N, ua, ub);
solveProgonka(delta, podDiag, mainDiag, nadDiag, funcF, N);
norma = norm1(delta,N);
//обновляем вектор решения
for (inti = 0; i<= N; i++)
{
z[i]=delta[i]+temp[i];
}
}
