Материал: Спецглавы высшей математики. численные методы. Пантелеев И.Н

Рис. IX.4. Пространственно-временная сетка, используемая для решения одномерного уравнения теплопроводности

Подставим данные выражения в исходное уравнение ut uxx и

разрешим получившееся уравнение относительно значений функции на верхнем временном слое. В результате получаем:

Формула (IX.10) решает поставленную задачу, поскольку она выражает решение в данный момент времени через решение в предыдущий момент времени.

Для решения смешенной краевой задачи необходимо аппроксимировать производную в граничном условии на правом конце:

ux (1, t) [u(1, t) g(t)].

Используя конечно-разностную аппроксимацию, получаем:

h1 [uin uin 1 ] [uin gi ].

208

Из (IX.11) находим:

Формулы (IX.10), (IX.12) позволяют начать вычисления по явной схеме.

Алгоритм вычислений по явной схеме реализуется следующей последовательностью действий:

1. Находим решение на сеточном слое t t , используя явную формулу:

поднявшись на одну строку вверх, т.е. увеличив i на единицу и используя uij с предыдущей строки. Аналогично вычисляется

решение в последующие моменты времени t 3 t,4 t,

У явной схемы имеется один существенный недостаток: если шаг по времени оказывается достаточно большим по сравнению с шагом по х, то погрешности округления могут стать настолько большими, что полученное решение теряет смысл, т. е. решение становится неустойчивым. Можно показать, что для применимости явной схемы должно выполняться условие k / h2 0.5.

Пример IX.2. Найти в пакете МАТLАВ решение краевой задачи уравнения теплопроводности (IX.8) с начальными условиями:

209

на временном интервале [0, 49].

Решение:

>>Nt=50; % число шагов по времени >>Nx=50; % Число узлов координатной сетки

>>t=1:Nt;

>>x=1:Nx;

% задание начальных и граничных условий

>>f(1,x)=0;

>>f(t,1)=0;

>>f(t,Nx)=0;

>>f(1,25)=1;

%вычисление решения в последовательные моменты времени в

%соотношении с (10.10)

>>for t=2:Nt

for x=2:Nx-1 f(t,x)=f(t-1,x)+0.15*(f(t-1,x-1)-2*f(t-1,x)+f(t-1,x+1));

end;

end;

>>surf(f) % визуализация численного решения (рис. IX.5)

Рис. IX.5. Численное решение уравнения теплопроводности

Получим явную вычислительную схему для нахождения численного решения УЧП гиперболического типа, на примере волнового уравнения:

Запишем уравнение (IX.15) в конечных разностях:

Полученное уравнение позволяет выразить значение функции и в момент времени ti 1 через значения функции в предыдущие моменты времени:

Разностная схема (IX.17) устойчива, если h / .

Отметим, что для начала вычислений в соответствии с (IX.17) необходимо знать значения функции в моменты времени t0 и t1 . Данное обстоятельство обусловлено тем, что мы решаем уравнение второго порядка по времени. Для нахождения решений в моменты времени t0 и t1 необходимо использовать, например, начальное условие для первых производных функции u (рис. IX.6).

Рис. IX.6. Решение волнового уравнения в моменты времени t=0, 9, 24, 49

Пример IX.3. Найти в пакете МАТLАВ решение краевой задачи волнового уравнения (IX.14), в котором 1,с начальными условиями:

на временном интервале [0,59].

Решение:

>>Nt=60; % число шагов по времени >>Nx=100; % число узлов по координате % задание начальных и граничных условий

>>U(1,j)=sin(2*pi*(j-1)/(Nx-1); % начальное условие

>>U(2,j)=U(1,j); % граничное условие

>>a=1;k=1;

>>for i=2:Nt

212

for j=2:Nx-1

U(i+1,j)=a.^2*k*(U(i,j+1)-2*U(i,j)+U(I,j-1))+2*U(I,j)-U(i-1,j); end;

end;

% визуализация решений уравнения в выбранные моменты времени

>>j=1:Nx; >>plot(2*pi*(j-1)/(Nx-1),U(1,j))

>>hold on

>>plot(2*pi*(j-1)/(Nx-1),U(10,j),':') >>plot(2*pi*(j-1)/(Nx-1),U(25,j),'--') >>plot(2*pi*(j-1)/(Nx-1),U(50,j),'-')

9.4. Неявная разностная схема для уравнения параболического типа

Построим неявную разностную схему для уравнения теплопроводности. Рассмотрим задачу нахождения решения УЧП

Воспользуемся следующими конечно-разностными аппроксимациями частных производных ut и uxx :

ut (x, t) k1 [u(x, t k) u(x, t)],

uxx (x, t) h 2 [u(x h, t k) 2u(x, t k) u(x h, t k)]

+ 1h 2 [u(x h, t) 2u(x, t) u(x h, t),

где λ – выбирается из отрезка [0,1].

213

Отметим, что uxx аппроксимируется взвешенным средним

центральных разностных производных в момент времени t и (t+k). При λ = 0.5 получаем обычное среднее двух центральных производных, при λ = 0 получаем обычную явную схему, о которой упоминалось ранее.

После замены частных производных ut , uxx в задаче (IX.18) получаем разностную задачу:

Перенесем все неизвестные значения и с верхнего временного слоя (с индексом (i+ 1)) в левую часть уравнения:

где r k / h 2 .

Схема уравнения системы (IX.20).

Если i фиксировано, а j изменяется от 2 до (n-1), соотношения (IX.20) определяют систему (n-2) уравнений с (n-2) неизвестными ui 1,2 , ui 1,3 , ui 1,4 , , ui 1,n 1 , которые являются решением задачи во внутренних узлах сетки на временном слое t=(i+1)∆t.

Наглядное представление о структуре каждого уравнения дает таблица.

где ui,1 = и(t,0), ui,n = и(t,1), которую будем использовать далее

при решении смешанной краевой задачи для уравнения параболического вида.

Пример IX.4. Найти в пакете МАTLАВ решение смешанной краевой задачи УЧП

и граничными условиями

u(t,0) 0, u(t,1) 1

с использованием неявной разности схемы (λ=1).

Для нахождения решения смешанной краевой задачи необходимо:

1. Создать файл U1.m (листинг IX.2), содержащий описание функции, возвращающей значения решения смешанной краевой задачи.

Листинг IX.2. Файл U1.m function z=U1(u,A,r,Nt,Nx)

%u - матрица, содержащая начальные и граничные условия

%А – матрица системы (IX.22)

%r=k/h2 – переменная, введенная в (IX.20)

%Nt – число шагов по времени

%Nx – число узлов по координате

z=u; u1=u(1:1,2:Nx-1); for i=1:Nt-1

a=r*z(i,1); b=r*z(i, Nx); u1(1,1)=u1(1,1)+a;

u1(1,Nx-2)=u1(1,Nx-2)+b; u2=A^-1*u1';

u1=u2';

for j=2:Nx-1 z(i+1,j)=u2(j-1,1);

end;

end;

217