Учебное пособие_Численные методы -26 -дек_2014 (97-2003)
.pdf
Решение в редакторе Visual Basic
Функции Urtem_1() вычисляет сеточную функцию ui,kи возвращает ее значения. Выделим диапазон ячеек B2:L22, введем формулу =Urtem_1() и нажмем комбинацию ctrl+shift+enter. Получим значения приближенного решения задачи. Вычислим в диапазоне B26:L46 значения точного решения. Для этого введем в строке B25:L25 значения переменной x: 0;0,1;…;1. В столбце A26:A46значения y:0;0,1;…;2. В ячейку B26 введем формулу =SIN(B$25)*EXP(-$A26) и протянем ячейку B26 маркером заполнения до L26, а затем протянем строку B26:L26 вниз до 46 строки. В ячейку N26 введем формулу =КОРЕНЬ(СУММКВРАЗН(B2:L22;B26:L46)) и вычислим сумму квадратов разностей значений. Получим значение
1,61293E-07.
Реализация в Visual Basic представлена на рис. 9.5 – рис. 9.6.
а)
Рис. 9.5. Решение задачи в Visual Basic: а) численно, б) в виде формул
271
б)
Рис. 9.5. (окончание)
а)
Рис. 9.6. Решение задачи в Visual Basic: а) численно, б) в виде формул
272
б)
Рис. 9.6. (окончание)
Код программы в редакторе Visual Basic:
Function mu(x): mu = Sin(x): End Function Function mu_1(t): mu_1 = 0: End Function
Function mu_2(t): mu_2 = Sin(1) * Exp(-t): End Function Function Urtem_1()
Application.Volatile (True)
Dim u(10, 20), fi(10), a(10), b(10) h = 0.1: tau = 0.1: N = 10: M = 20
For i = 0 To 10: u(i, 0) = mu(i * h): Next i For k = 0 To M - 1
tk = tau * k: u(0, k + 1) = mu_1(tk + tau): u(N, k + 1) = mu_2(tk + tau) For i = 1 To N - 1
fi(i) = (2 - 2 * h ^ 2 / tau) * u(i, k) - u(i - 1, k) - u(i + 1, k) Next i
d = 2 + 2 * h ^ 2 / tau: a(1) = 1 / d: b(1) = -(fi(1) - mu_1(tk + tau)) / d For i = 2 To N - 2
a(i) = -1 / (a(i - 1) - d): b(i) = (fi(i) - b(i - 1)) / (a(i - 1) - d) Next i
u(N - 1, k + 1) = (fi(N - 1) - mu_2(tk + tau) - b(N - 2)) / (a(N - 2) - d) For i = N - 2 To 1 Step -1
u(i, k + 1) = a(i) * u(i + 1, k + 1) + b(i) Next i
Next k
Urtem_1 = Application.Transpose(u) EndFunction
273
Пример 9.4
Найти значения функции u(x,t) , удовлетворяющие уравнению ut |
uxx |
с |
начальными условиями (т.е. при t=0) |
u(x,0) sin(πx), 0 x 1 |
, краевыми |
условиями первого рода u(0,t) 0, u(1,t) 0 . |
|
|
Выберем по аргументу x шаг h=0,1, т.е. n=10. Так как q=0,5, то по аргументу t получаем шаг k=h2/2=0,005.
Решение в Microsoft Excel.
Откроем новый рабочий лист Excel. Отведем столбец А под номера временных слоев от 0 до, например, 25. Заполним блок А3:А28 числами от 0 до 25. Во второй строке в ячейки В2:L2 занесем значения пространственной координаты Х от 0 до 1. В блок В3:В28 и в блок L3:L28 занесем нули - значения краевых условий. В ячейку С3 занесем формулу =SIN(3,1415*C2) и скопируем ее в блок D3:K3. Тем самым введем в нулевой временной слой начальные условия. Далее в ячейку С4 введем основную расчетную формулу, полученную из алгебраического конечно-разностного уравнения =(B3+D3)/2. Эту формулу надо скопировать в блок С4:L28.
Результаты вычислений представлены в таблице (рис. 9.7) для первых 19 слоев.
Таблица не наглядна, а при увеличении количества временных слоев вообще становится нечитаемой. Результаты расчетов лучше всего представить в виде диаграммы. Для построения диаграммы надо выделить блок А2:L28 и вызвать Мастер Диаграмм. Появится диалог из пяти шагов. На шаге 1 из 5 следует убедиться в том, что выделен блок А2:L28 и щелкнуть по кнопке Шаг>. На шаге 2 из 5 нужно выбрать тип диаграммы - Поверхность. На шаге 3 из 5 выберем формат 1 диаграммы. Далее на шаге 4 из 5 выберем Ряды данных - в строках и в полях
Отвести 0 строк введем 1 так, чтобы получилось Отвести 1 строку для меток оси Х.
274
Рис. 9.7. Результаты вычислений
То же самое проделаем для оси Y - Отвести 1 столбцов для меток оси Y. На последнем шаге 5 из 5 в поле Название диаграммы введем текст «Уравнение теплопроводности», в поле Название по оси категорий (Х) введем текст – «длина», в поле Значений (Z) текст «температура», в поле Рядов (Y) текст «слой». Щелкнув по кнопке Закончить, получим диаграмму, примерный вид которой приведен на рис. 9.8.
275
Р9u(i, j, 0) = mu(x, y): Next j: Next i '1.k>0
For k = 0 To M - 1 '1
tk = tau * k: tk1 = tk + tau: tk12 = tk + tau / 2 For j = 1 To Ny - 1: y = hy * j
'1.1
For i = 1 To Nx - 1: x = hx * i
fi(i) = hxy * u(i, j, k) - hx2 * ((u(i, j - 1, k) + u(i, j + 1, k)) / hy2 + f(x, y, tk12)) Next i
'1.2
d = 2 + 2 * hx2 / tau
m10 = (mu_1(y, tk1) + mu_1(y, tk)) / 2
m10 = m10 - (tau / 4) * (mu_1(y - hy, tk1) - 2 * mu_1(y, tk1) + mu_1(y + hy, tk1)) / hy2 m10 = m10 + (tau / 4) * (mu_1(y - hy, tk) - 2 * mu_1(y, tk) + mu_1(y + hy, tk)) / hy2 m2N = (mu_2(y, tk1) + mu_2(y, tk)) / 2
m2N = m2N - (tau / 4) * (mu_2(y - hy, tk1) - 2 * mu_2(y, tk1) + mu_2(y + hy, tk1)) / hy2 m2N = m2N + (tau / 4) * (mu_2(y - hy, tk) - 2 * mu_2(y, tk) + mu_2(y + hy, tk)) / hy2 a(1) = 1 / d: b(1) = (fi(1) - m10) / (-d)
For i = 2 To Nx - 2 a(i) = -1 / (a(i - 1) - d):
b(i) = (fi(i) - b(i - 1)) / (a(i - 1) - d) Next i
'1.3
uk12(Nx - 1, j) = (fi(Nx - 1) - m2N - b(Nx - 2)) / (a(Nx - 2) - d) For i = Nx - 2 To 1 Step -1: uk12(i, j) = a(i) * uk12(i + 1, j) + b(i): Next i
uk12(0, j) = m10: uk12(Nx, j) = m2N Next j
'2
For i = 1 To Nx - 1: x = hx * i '2.1
For j = 1 To Ny - 1: y = j * hy
fi(j) = hyx * uk12(i, j) - hy2 * ((uk12(i - 1, j) + uk12(i + 1, j)) / hx2 + f(x, y, tk12)) Next j
'2.2
d = 2 + 2 * hy2 / tau
a(1) = 1 / d: b(1) = (fi(1) - mu_3(x, tk1)) / (-d) For j = 2 To Ny - 2
a(j) = -1 / (a(j - 1) - d): b(j) = (fi(j) - b(j - 1)) / (a(j - 1) - d) Next j
'2.3
u(i, Ny - 1, k + 1) = (fi(Ny - 1) - mu_4(x, tk1) - b(Ny - 2)) / (a(Ny - 2) - d)
For j = Ny - 2 To 1 Step -1: u(i, j, k + 1) = a(j) * u(i, j + 1, k + 1) + b(j): Next j Next i
For j = 0 To Ny: y = hy * j
u(0, j, k + 1) = mu_1(y, tk1): u(Nx, j, k + 1) = mu_2(y, tk1) Next j
For i = 0 To Nx: x = hx * i
u(i, 0, k + 1) = mu_3(x, tk1): u(i, Ny, k + 1) = mu_4(x, tk1) Next i
278
Next k
For i = 0 To 10: For j = 0 To 10: uy(i, j) = u(i, j, M): Next j: Next i
Urtem_2 = Application.Transpose(uy) End Function
Полученное решение представлено на рис. 9.9.
Рис. 9.9. Метод сеток для решения эллиптических уравнений
279
Контрольные вопросы
1.Что такое дифференциальное уравнение в частных производных?
2.Перечислите классы уравнений в частных производных.
3.Какие задачи относятся к решению уравнений гиперболического, эллиптического, параболического типа?
4.Перечислите виды граничных условий для уравнений в частных производных.
5.Как можно задать начальное условие для дифференциального уравнения в частных производных?
6.Назовите аналитические методы решения дифференциальных уравнений в частных производных.
7.В чем отличие явной от неявной разностной схемы?
8.В чем суть метода сеток для решения параболических уравнений в частных производных?
9.В чем суть метода сеток для решения гиперболических уравнений?
10.В чем суть метода сеток для решения эллиптических уравнений?
280