Порядок выполнения лабораторной работы. Рецензент: А.А. Викарчук – профессор, заведующий кафедрой «Нанотехнологий и новых материалов», д

РЕШЕНИЕ КРАЕВОЙ ЗАДАЧИ ДЛЯ УРАВНЕНИЯ ПАРАБОЛИЧЕСКОГО ТИПА РАЗНОСТНЫМ МЕТОДОМ

Тольятти 2010

УДК

ББК

М

Рецензент: А.А. Викарчук – профессор, заведующий кафедрой «Нанотехнологий и новых материалов», д. ф.-м.н.

Методические указания к лабораторной работе № 5 «Решение краевой задачи для уравнения параболического типа разностным методом», Сафронов А.И., 2010 – 10с.

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

УДК

ББК

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

©Авторы пособия

Цель работы: усвоить численный метод решения начально - краевой задачи для уравнения параболического типа

Приборы и принадлежности: миллиметровая бумага.

ТЕОРЕТИЧЕСКИЙ МАТЕРИАЛ

Рассмотрим смешанную задачу для однородного уравнения теплопроводности. Задача состоит в отыскании функции и (х, t), удовлетворяющей в областиD ={ (х, t)|0 < х < а,0 <t < Т}уравнению (1) начальному условию u(x,t) = f(x) (2) и краевым условиям третьего рода

. (3)

К задаче (1)—(3) приводит, а частности, задача о распространении тепла в однородном стержне длины а, на левом конце которого задаётся тепловой поток в стержнь (q), а правый конец теплоизолирован.

В данной лабораторной работе рассматриваются краевые условия второго рода. Замена переменных t → kt приводит уравнение (1) к виду

поэтому вдальнейшем будем считать k=1.

Построимв области D = { (х, t) )|0 < х < а, 0 <t < Т}равномерную прямоугольную сетку а шагом h в направлении х и шагом τ — в направлении t (рис. 1). Обозначимузлы cетки через (x₁,t₁), …(x_i,t_j), …, а приближенные значения функции u(x, t) в узлах – u_ij. Тогда x_i=ih, i=0,1,…,n, h=a/n, t_j= jτ, j=0,1,…,m, τ=T/m.

Аппроксимируем уравнение (1) на четырехточечном шаблоне, который изображен на рис. 1 жирными линиями. В результате получаем неявную двухслойную разностную схему:

λu_i+1,_j — (1+2λ)u_{i, j} + λu_i-1,_j = – u_{i,j –1} (4)

которая аппроксимирует уравнение (1) с погрешнос тью О(τ + h²). Здесь λ = τ / h².

Рис. 1 Рис. 2

Схема (4) аппроксимирует уравнение (1) только во внутренних узлах сетки, поэтому число уравнений в схеме (4) меньше числа неизвестных u_ij. Недостающие уравнения получают из краевых условий

. ( 5 )

Схема (4)—(5) неявная, поэтому значения u_ij находят как решение системы линейных уравнений (4). Для решения системы (4) можно применять любой алгоритм решения систем линейных уравнений, однако система (4) обладает трехдиагональной матрицей и рациональнее всего решать ее методом прогонки. Таким образом, решив систему разностных уравнений, найдем значения функции и на временном слое j, если известно решение на временном слое j – 1.

Итак, алгоритм численного решения задачи имеет следующий вид. На нулевом временном слое (j = 0) решение известно из начального условия u_i0 = f (x_i). На каждом следующем слое искомая функция определяется как решение системы (4)—(5).

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

Уравнение (4) в общем виде записывается

a_i u_i-1 — 2b_i u_i + c_i u_i+1 = f_i , i= 1, 2, … , n-1 (6)

на левой границе области

c₀u₁ – 2b₀u₀ = f₀ (7)

на правой границе области

a_Nu_N-1-2b_Nu_N=f_N. (8)

Имеем u₀, u₁ , u₂ ,...

Будем искать решение в виде:

u_i-1=Р_iu_i+Q_i( * )

u₀=c₀/2b₀·u₁+(-f₀/2b₀) (9)

Р₁ = c₀ /2b₀ , Q₁ = - f₀ / 2b₀ ;

Для всех остальных Р_i , Q_i получим рекуррентные формулы, подставив (*) в (6):

a_i (Р_i u_i + Q_i) — 2b_i u_i + c_i u_i+1 = f_i ,

(a_i Р_i –2b_i ) u_i + c_i u_i+1 = f_i – a_i Q_i

u_i =

u_i = Р_i+1u_i+1 + Q_i+1 ,

Если известны Р₁, Q₁, то можно найти Р₂, Q₂, и т.д.

u_N-1 = P_N u_N +Q _N ;

u_N-1 = 2b_N /a_N ·u_N + f_N /a_N ;

Имеем систему двух уравнений.

Отсюда: u_N =

Зная u_N по формуле (*) вычисляем все u_iот i = N-1 до 0.

Таким образом, находятся все u_i на верхнем слое.

Замечательным свойством неявной схемы (4) является ее устойчи­вость при любых значениях параметра λ = τ / h² > 0. Преимущества схемы (4) особенно ощутимы при сравнении с явной схемой, которая
получается при аппроксимации уравнения (1) на шаблоне, изображен­ном на рис. 2.

Явная схема оказывается устойчивой только при λ < 1/2, т. е. при τ < h² /2. Это означает, что вычисления по явной схеме придется вести с очень малым шагом по τ, что конечно может привести к большим за­тратам машинного времени. В неявной схеме вычисления на одном шаге требуют больше операций, чем в явной схеме, но зато величину шага τ можно выбирать, как угодно большой без риска нарушить устойчивость схемы. Все это позволяет значительно уменьшить машинное время, необходимое для решения задачи. Схема (4) обладает сходимостью. Это означает, что при h, τ → 0 решение разностной задачи (4) — (5) стремится к точному решению смешанной задачи (1)—(3).

Для численного решения задачи (1)—(3) предназначе на подпрог­рамма CALC.

Подпрограмма CALC решает систему (4) методом прогонки и вычисляет решение на каждом слое по значениям решения с предыдущего слоя.

Входные параметры: НХ — шаг сетки h по переменной х; НТ — шаг сетки τего переменной t; N — количество шагов сетки по х; U — массив из N + 1 чисел, содержащий перед обращением к CALC начальное значение решения в узлах первого слоя сетки; Р, Q — рабочие массивы из N — 1 чисел каждый /это прогоночные коэффициенты/.

Выходные параметры: U — массив из N + 1 чисел (он же входной), содержащий значение решения в узлах очередного временного слоя сетки.

Перед обращением к программе CALC необходимо:

1)описать массивы U, P, Q:

2)присвоить значения параметрам НХ, НТ, N;

3)присвоить элементам U начальное значение решения в узлах первого слоя сетки.

Задание. Используя программу CALC, решить краевую зада­чу с граничными условиями 2-го рода для однородного уравнения теплопроводности в области D = { (х, t)|0 < х < а, 0 <t < Т}.Результаты печатать на каждом временном слое.

Порядок выполнения лабораторной работы

1.Составить основную программу, содержащую цикличес кое обращение к подпрограмме решения параболического уравнения CALC с печатью результатов на каждом шаге. Образец дан в приложении1.

2. Составить подпрограмму CALC решения параболичес кого уравнения методом прогонки.

3. Составить подпрограмму - функцию, вычисляющую значение функции F(x) для начального значения u( x, 0 ).

4. Провести вычисления на компьютере.

5. Результаты решения нанести на мимиллиметровку в виде графика изменения функции u по времени.

Примечание: для более плавного поведения функции увеличить начальное количество точек разбиения отрезка X [ 0 , 1 ] в два раза.

Варианты заданий

Решить краевую задачу для уравнения теплопровод ности с начальным условием u(х, 0) = f (х)и граничными условиями . Во всех вариантах q=1000, Ψ=100, α = 100.

РЕКОМЕНДУЕМАЯ ЛИТЕРАТУРА

1. А.А.Самарский, А.В.Гулин. Численные методы. М.: Главная ред. Физ.- мат. лит. - 1989 г., 430 с.

2. Бахвалов Н.С. Численные методы / Н. С.Бахвалов, Н. П. Жидков, Г.М. Кобельков. – изд.4-е; Гриф МО. – М.: Бином. Лаборатория знаний. 2006 г. – 636 с.

Приложение1

Листинг основного блока программы на Pascal

program PAR_5;

const ht=0.02;

n=10;

hx=1/n;

var j,i,i1,N1,N2:integer;

x,y,t,T1,al,B1,A1:real;

u:array [0..n] of real;

P:array [0..n] of real;

Q:array [0..n] of real;

XP:array [0..n] of real;

begin

XP[0]:=0.0;

XP[10]:=1.0;

u[0]:=10;

u[10]:=20;

T1:=9*HT;

T:=0;

for i:=1 to n do

begin

x:=i*HX;XP[i]:=x;

u[i]:=f(x);

end;

writeln(T:8:5);

for i:=0 to n do write(u[i]:6:3);

writeln;

for j:=1 to 40 do

begin

calc; T:=T+HT;

if (T>T1) then begin writeln(T:8:5);

for i:=0 to n do write(u[i]:6:3);

writeln;

T1:=T1+10*HT;end;

end;

end.

Методические указания

по выполнению лабораторной работы № 5