Задачи №2 и 13

Первое задание

Задача №2

Постановка задачи:

Используя схему бегущего счета и итерационные методы решить задачу:

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

Построим характеристики для t₀=0 и для различных x₀ (сплошные линии) и для x₀=0 и различных t₀ (пунктирные линии). При этом уравнение примет вид

В результате получим следующее семейство характеристик :

Как видно, на интервале характеристики не пересекаются, следовательно нет опрокидывания волны. Решение будем строить на временном интервале от 0 до 2.

Перейдем теперь к построению разностной схемы. Введем разностную сетку:

Здесь N – число узлов по х, τ – шаг по времени.

Перепишем первое уравнение исходной системы в дивергентном виде:

Введем сеточную функцию y_ij=ũ(x_i,y_j). Разностная аппроксимация уравнения в точке (x_i+h_x/2,t_j+τ/2) выглядит следующим образом:

Для граничных и начальных условий:

В данном случае используется четырехточечный шаблон:

Раскладывая u(x,t) в ряд Тейлора в центре ячейки можно получить, что невязка равна O(τ²+h²), то есть порядок аппроксимации данной схемы равен 2.

Устойчивость схемы доказывается с помощью принципа максимума. Запишем двуслойную разностную схему в виде

, где суммирование на каждом слое проводится по узлам шаблона около n-го узла. Коэффициенты перенумеруем так, что . Тогда схема равномерно устойчива по начальным данным, если

, C = const

И устойчива по правой части, если

Зафиксируем решение на исходном слое и внесём погрешность в правую часть. Тогда погрешность на новом слое

Отсюда

выбираем узел и заменяем все величины их максимумами

Отсюда с учётом второго неравенства следует, что , то есть схема устойчива по правой части.

В рассматриваемой задаче отсутствует неоднородность, поэтому нет необходимости проверять устойчивость по правой части. Таким образом, достаточно проверить устойчивость по начальным данным. Переписывая нашу схему в форме двух сумм получаем

при

при n = 0 и n = N

Видно, что условие устойчивости выполнено во всех узлах сетки, то есть схема безусловно устойчива по начальным данным и краевым условиям.

Из аппроксимации и устойчивости следует сходимость разностной схемы.

Введем функцию ):

Значение функции ŷ_i+1 находится методом последовательных приближений, при этом первым приближением считается значение ŷ_i. Значение же ŷ_i+1 находится по формуле

Решение строилось на интервалах:

Подсчеты велись до достижения точности 1%, либо до 1000 итераций.

В результате был получен следующий график:

Текст программы (MathLab):

function OMM1

%Основной m-файл

%Количество шагов

nx=100;

nt=200;

%Ограничение по времени

tmax=2;

%Шаг х

dx=-1/nx;

%Шаг t

dt=tmax/nt;

%Точность для решения уравнения на "y" методом касательных

precision=0.01;

%Для графика: оси x и t с делениями

x=0:dx:-1;

t=0:dt:tmax;

%Граничное условие

for k=0:1:nx

y(k+1,1)=2-4*atan(dx*k+2)/3.141592636;

gr(k+1)=y(k+1,1);

end;

%Начальное условие

for k=0:1:nt

y(1,k+1)=(2-4*atan(2)/3.141592636)*exp(-(dt*k));

nach(k+1)=y(1,k+1);

end;

%Поиск y(i+1, j+1):

for m=1:nt

for n=1:nx

y(n+1,m+1)=res(y(n,m),y(n+1,m),y(n,m+1),dx,dt,precision);

end;

end;

surf(t,x,y);

end

function u=res(u1,u2,u3,dx,dt,eps)

%Это чтобы хотябы один шаг был выполнен

deltaUp=eps+1;

%Начальное приближение

aproxy2=u2;

l=1;

while deltaUp>eps

l=l+1;

%aproxy2 - это y(i+1, j+1, s+1), a aproxy - это y(i+1, j+1, s)

%(см.описание)

aproxy=aproxy2;

%Схема метода касательных: y(s+1)=y(s)-f(y(s))/f'(y(s))

%u1 - это y(i,j)

%u2 - это y(i+1,j)

%u3 - это y(i,j+1)

f=-aproxy^2/(4*dx)+aproxy/(2*dt)+(u3-u1-u2)/(2*dt)-(u2^2-u1^2-u3^2)/(4*dx);

prf=-aproxy/(2*dx)+1/(2*dt);

aproxy2=aproxy-f/prf;

deltaU=abs(aproxy2-aproxy);

deltaUp = deltaU;

if(l>10000)

deltaUp = 0.001;

r = deltaU;

end

end;

u=aproxy2;

end

Второе задание

Задача №13

Постановка задачи:

Решить краевую задачу методом переменных направлений.

Первым шагом вводится сетку в области :

Вторым шагом необходимо ввести разностную аппроксимацию оператора Лапласа.

При этом из значения функции на нулевом слое известно из начальных условий.

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

То есть мы вводим промежуточный слой k+1/2, в результате переход со слоя k на слой k+1 будет осуществляться в два этапа: сначала решается первое уравнение, неявное по направлению x и явное по направлению y, а затем второе, явное по направлению x и неявное по направлению y.

Невязка для данной схемы равна O(|h|2+τ2), то есть ее порядок аппроксимации равен 2. Кроме того, схема безусловно устойчива (см. Самарский А. А. «Теория разностных схем»). Следовательно, схема сходится.

Используя явный вид разностных операторов, приходим к задаче:

Эта задача состоит из двух подзадач: переход со слоя k на слой k+1/2 и переход со слоя k+1/2 на слой k+1. Каждая из этих задач решается методом прогонки. Первая задача решается для всех j=1…N2‑1, а вторая – для всех i=1…N1‑1. Рассмотрим этот метод поподробнее. Обе наши задачи имеют следующий общий вид:

_,

_{причем из граничных условий получаем}

Будем искать решение этой задачи в виде ωi=Ai+1ωi+1+Bi+1. Подставляя данный вид в общий вид задачи и исключая ωi, получим рекуррентное соотношение для коэффициентов прогонки:

Вычисление ωn:

Выразим ωn-1 через ωn

в результате подстановки получим:

Далее задача решается в два этапа:

• Сначала определяются все коэффициенты прогонки, счет идет от 2 до N1,2,

• Зная коэффициенты прогонки определяем значения функции в узлах, счет идет от N1,2-1 до 1.

Метод прогонки устойчив при |Сi|≥|Ai|+|Bi| , в нашем случае всегда выполняется строгое равенство.

В результате были получен следующий набор графиков:

t = 1 t=2

y =0.2

y=0.5

x =0.3

x=0.5

Текст программы (MathLab):

function nomm2

global pi;

pi=3.141592636;

% node

N1=101;

N2=101;

J=201;

% initial intervals and steps

T=2;

h1=1/(N1-1);

h2=1/(N2-1);

tau=T/(J-1);

x=0:h1:1;

y=0:h2:1;

t=0:tau:T;

% initial time layer

for i=1:N1;

for j=1:N2;

u(1,i,j)=0;

end

end

for k=2:J;

% boundary conditions

for i=1:N1;

w(i,1)=0;

w(i,N2)=0;

u(k,i,1)=0;

u(k,i,N2)=0;

end

for j = 1:N2;

w(1,j)=0;

w(N1,j)=0;

u(k,1,j)=0;

u(k,N1,j)=0;

end

% layer k-1+tau/2

for j=2:N2-1;

for i=2:N1-1;

a(i-1)=0.5*tau/(h1^2);

b(i-1)=-(1+tau/(h1^2));

c(i-1)=a(i-1);

d(i-1)=-(0.5*tau/(h2^2)*(u(k-1,i, j-1)+u(k-1,i,j+1))+
(1-tau/(h2^2))*u(k-1,i,j)+0.5*tau*sin(pi*h1*(i-1))
*h2*(j-1)*(h2*(j-1)-1)*tau*(k-1));

end

v=sweep(a,b,c,d,N1-2);

for i=2:N1-1;

w(i,j)=v(i-1);

end

end

% layer k

for i=2:N1-1;

for j=2:N2-1;

aa(j-1)=0.5*tau/(h2^2);

bb(j-1)=-(1+tau/(h2^2));

cc(j-1)=aa(j-1);

dd(j-1)=-(0.5*tau/(h1^2)*(w(i-1,j)+w(i+1,j))+(1-tau/
(h1^2))*w(i,j)+0.5*tau*sin(pi*h1*(i-1))*h2*
(j-1)*(h2*(j-1)-1)*tau*(k-0.5));

end

vv=sweep(aa,bb,cc,dd,N2-2);

for j=2:N2-1;

u(k,i,j)=vv(j-1);

end

end

end

% Вывод результатов

k = 101;

for j=1:N2;

for i=1:N1;

v(i,j)=u(k,i,j);

end

end

surf(x,y,v);

pause;

k = 201;

for j=1:N2;

for i=1:N1;

v(i,j)=u(k,i,j);

end

end

surf(x,y,v);

pause;

j=21;

for k=1:J;

for i=1:N1;

f(k,i)=u(k,i,j);

end

end

surf(x,t,f);

pause;

j=51;

for k=1:J;

for i=1:N1;

f(k,i)=u(k,i,j);

end

end

surf(x,t,f);

pause;

i=31;

for j=1:N2;

for k=1:J;

r(k,j)=u(k,i,j);

end

end

surf(y,t,r);

pause;

i=51;

for j=1:N2;

for k=1:J;

r(k,j)=u(k,i,j);

end

end

surf(y,t,r);

end

function sweep = sweep(a, b, c, d, N)

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

A(1)=-c(1)/b(1);

B(1)=d(1)/b(1);

% direct sweep

for i=2:N-1;

A(i)=-c(i)/(a(i)*A(i-1)+b(i));

B(i)=(d(i)-a(i)*B(i-1))/(a(i)*A(i-1)+b(i));

end

x(N)=(d(N)-a(N)*B(N-1))/(b(N)+a(N)*A(N-1));

% back sweep

for i=N-1:(-1):1;

x(i)=A(i)*x(i+1)+B(i);

end

sweep=x;

end