Решение задачи. Классификация краевой задачи и её физический смысл

Классификация краевой задачи и её физический смысл

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

(1)

Произведём замену переменных: . В новых переменных уравнение приобретает вид:

Штрихи в дальнейшем будем опускать.

Линейное дифференциальное уравнение второго порядка такого вида относится к уравнениям гиперболического типа. Данное уравнение описывает колебание струны в зависимости от времени. Здесь переменная t имеет физический смысл времени.

Первая краевая задача для уравнения колебания струны состоит в поиске решения уравнения, удовлетворяющего начальному и граничным условиям.

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

Начальные условия (в начальный момент времени t=0):

- задаёт форму струны.

- задаёт распределение скоростей.

Граничные условия:

- означает, что по такому закону движется левый конец струны.

- вид граничного условия третьего рода на правом конце струны.

-плотность внешней силы, рассчитанная на единицу длины.

РЕШЕНИЕ ЗАДАЧИ С ПОМОЩЬЮ ЯВНОЙ РАЗНОСТНОЙ СХЕМЫ

j+1 i i+1

j

j-1

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

Конечно-разностная схема для данной задачи имеет вид:

с погрешностью аппроксимации порядка .

Обозначим и выразим через остальные значения сеточной функции, входящие в уравнение:

(3)

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

Аппроксимация 1-го начального и 1-го граничного условий:

(4)

(5)

Аппроксимация 2-го начального условия:

Для аппроксимации 2-го начального условия (с погрешностью порядка ) разложим функцию U(x,t) в окрестности точки (x,0) по формуле Тейлора:

.

Учитывая уравнение краевой задачи (1) и условие (2) запишем:

.

Перейдём к конечным разностям, записываемым в узле (i,1) (т.е. на первом временном слое):

(6)

Данная формула используется на начальном этапе для вычисления значения функции U(x,t) на первом слое, по известным значениям функции на нулевом слое и на границе.

Аппроксимация 2-го граничного условия:

Для аппроксимирования 2-го граничного условия (с погрешностью порядка ) разложим функцию U(x,t) в окрестности точки ( 1 ,t) по формуле Тейлора:

Используя уравнение краевой задачи (1) и условие (2) получаем:

Перейдя к конечным разностям, записываемым в узле (M-1,j), получаем:

(7)

Выразим из (7) :

(8)

Подставляя в (7) j=0 и учитывая (4), получаем:

(9)

Подставляя далее в (6) i=M и учитывая (4), получаем:

(10)

Из (3) при i=M, j=0, получим:

Учитывая (4) и подставляя сюда (9) и (10), находим:

(11)

Объединяя полученные результаты, получаем аппроксимацию краевого условия третьего рода, т.е. значение функции на границе i=M, в явном виде:

Текст функции, вычисляющий матрицу решения размерностью MxN, приведены в Приложении.

Устойчивость решения:

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

,

т.е. если это условие устойчивости не будет выполнено, то в процессе рекуррентного решения возможно накапливание ошибок от слоя к слою. Отсюда, в частности, получаем для явной схемы () условие устойчивости Куранта-Леви:

. (проверка этого условия реализована в программе)

Итак, U(i,j) при j=0 и j=1 определены. Включается рекуррентная процедура описываемая уравнением (3) и вычисляется U(i,j) для всех i=1,2,...,M-1, для каждого фиксированного j=2,...,N-1.

РЕШЕНИЕ ЗАДАЧИ С ПОМОЩЬЮ НЕЯВНОЙ РАЗНОСТНОЙ СХЕМЫ

Для аппроксимации используем Т-образный пятиточечный шаблон:

Уравнение аппроксимируется следующей системой разностных уравнений:

(12)

с погрешностью аппроксимации .

Обозначим и запишем (12) следующим образом:

Дополним это уравнение уравнениями аппроксимации начальных и граничных условий.

Аппроксимация 1-го начального и 1-го граничного условий:

Аппроксимация 2-го начального условия:

Для аппроксимирования 2-го начального условия (с погрешностью порядка ) разложим функцию U(x,t) в окрестности точки (x,0) по формуле Тейлора:

.

Учитывая уравнение краевой задачи (1) и условие (2) и переходя к конечным разностям:

(13)

Эта формула отличается от аналогичной для явной схемы тем, что аппроксимация разностной производной второго порядка по x производится на первом слое, а не на нулевом. Запишем (13) к виду удобному для применения метода прогонки:

(14)

Аппроксимация 2-го граничного условия:

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

Преобразуем их к виду, удобному для использования метода прогонки:

Вычисление прогоночных коэффициентов.

Сначала найдем U(i,j) на слое j=1. Определим прогоночные коэффициенты. Учитывая 1-ое граничное условие и уравнение (14) получаем:

Теперь вычислим граничные прогоночные коэффициенты на том же слое:

Методом прогонки находим U(i,1) где i=1...M;

Теперь зная значения функции при j=0,1, мы можем найти U(i,j), где j=2...N.

Используя уравнение (2) и 1-ое граничное условие находим прогоночные коэффициенты:

Теперь вычислим граничные прогоночные коэффициенты:

Методом прогонки находим U(i,j), где i=1...M, j=2..N.

СРАВНЕНИЕ ЯВНОЙ И НЕЯВНОЙ СХЕМ:

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

Рассмотрим матрицы следующих размерностей: 250х250, 500х500, 1000х1000, 2000х2000.

Таблица норм сравнения явной и неявной схем.

Размерность сетки	2000х2000	1000х1000	500х500	250х250
Макс. расхождение	0.1245	0.2498	0.5029	1.0190

Решение модельной задачи:

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

решением которой является функция .

Сравнение графиков для явной и неявной схем:

Результаты решения задачи методом явной схемы на сетке 1000х1000:

Результаты решения задачи методом неявной схемы на сетке 1000х1000:

Результаты решения задачи методом явной схемы на сетке 1000х1000:

Результаты решения задачи методом неявной схемы на сетке 1000х1000:

Результаты решения задачи методом явной схемы на сетке 1000х1000:

Результаты решения задачи методом неявной схемы на сетке 1000х1000:

Результаты решения модельной задачи на сетке 1000х1000:

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

Размерность сетки	2000х2000	1000х1000	500х500	250х250
Макс. ошибка	4.814089e-007	1.925694e-006	7.700859e-006	3.078770e-005

Таким образом, при увеличении разбиения сетки ошибка численного решения уменьшается.

Матрица решений полученная явным методом разносного решения на сетке 1000х1000:

x=0	271.8	495.3	902.5	1644.5	2996.4	5459.8
x=0.1	222.6	405.5	738.9	1346.4	2453.3	4470.1
x=0.2	182.2	332.0	605.0	1102.3	2008.6	3659.8
x=0.3	149.2	271.8	495.3	902.5	1644.5	2996.4
x=0.4	122.1	222.6	405.5	738.9	1346.4	2453.3
x=0.5	100.0	182.2	332.0	605.0	1102.3	2008.6
x=0.6	81.9	149.2	271.8	495.3	902.5	1644.5
x=0.7	67.0	122.1	222.6	405.5	738.9	1346.4
x=0.8	54.9	100.0	182.2	332.0	605.0	1102.3
x=0.9	44.9	81.9	149.2	271.8	495.3	902.5
x=1.0	36.8	67.0	122.1	222.6	405.5	738.9

Матрица решений полученная неявным методом разносного решения на сетке 1000х1000:

x=0	271.8	495.3	902.5	1644.5	2996.4	5459.8
x=0.1	222.6	405.5	738.9	1346.4	2453.4	4470.3
x=0.2	182.2	332.0	605.0	1102.4	2008.7	3660.1
x=0.3	149.2	271.8	495.3	902.5	1644.6	2996.6
x=0.4	122.1	222.6	405.5	738.9	1346.5	2453.4
x=0.5	100.0	182.2	332.0	605.0	1102.4	2008.7
x=0.6	81.9	149.2	271.8	495.3	902.6	1644.6
x=0.7	67.0	122.1	222.6	405.5	738.9	1346.5
x=0.8	54.9	100.0	182.2	332.0	605.0	1102.4
x=0.9	44.9	81.9	149.2	271.8	495.3	902.6
x=1.0	36.8	67.0	122.1	222.6	405.5	739.0

Текст функции, реализующий решение задачи методом явной схемы:

clear;

t=1;

h=0;

while (t>h)

M = input('Введите число узлов по оси X:');

N = input('Введите число узлов по оси T:');

A=1;

B=1/4;

h=A/M; %длина шага разбиения по оси X

t=B/N; %длина шага разбиения по оси Т

if (t>h) disp('НЕ ВЫПОЛНЕНО УСЛОВИЕ СХОДИМОСТИ'); end;

end;

g=(t/h)^2;

U=zeros(M+1,N+1); %матрица решений

f=zeros(M+1,N+1); %матрица значений функции fi,j

x = zeros(M+1,1);

y = zeros(N+1,1);

for i = 1:(M+1),

x(i) = (i-1)*h;

end;

for j = 1:(N+1),

y(j) = (j-1)*t;

end;

p=zeros(M+1,1); %первое начальное условие

q=zeros(M+1,1); %второе начальное условие

fi=zeros(N+1,1); %первое граничное условие

for i=1:M+1

p(i)=100*exp(1-2*x(i));

q(i)=1200*exp(1-2*x(i));

U(i,1)=p(i);

end;

for j=1:N+1

fi(j)=100*exp(12*y(j)+1);

U(1,j)=fi(j);

end;

%формирование матрицы значений функции f(i,j)

for i=1:M+1

for j=1:N+1

f(i,j)=14000*exp(12*y(j)-2*x(i)+1);

end;

end;

%заполнение первого слоя по значениям на нулевом

for i=2:M

U(i,2)=p(i)+t*q(i)+((t^2)/2)*((p(i-1)-2*p(i)+p(i+1))/(h^2)+f(i,1));

end;

U(M+1,2)=g*(p(M)-p(M+1)+h*(-2*U(M+1,1)))+p(M+1)+t*q(M+1)+(t^2/2)*f(M+1,1);

%основной цикл – вычисление массива решений

for j=2:N

for i=2:M

U(i,j+1)=2*(1-g)*U(i,j)+g*(U(i-1,j)+U(i+1,j))-U(i,j-1)+t^2*f(i,j);

end;

U(M+1,j+1)=2*g*(U(M,j)-U(M+1,j)+h*(-2*U(M+1,j)))+f(M+1,j)*t^2+2*U(M+1,j)-U(M+1,j-1);

end;

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

y=zeros(M+1,11);

for i=1:M+1

for j=1:11

y(i,j)=U(i,(j-1)*N/10+1);

end;

end;

i=0:1:M;

plot(i/M,y);

legend('t=0','t=1/10 T','t=2/10 T','t=3/10 T','t=4/10 T','t=5/10 T','t=6/10 T','t=7/10 T','t=8/10 T','t=9/10 T','t=T','Location','NorthEastOutside')

title('Параметрический график при различных t');

xlabel('Ось X ');

ylabel('Ось U ');

pause;

X = 0:h:A; %массив делений по оси X

Y = 0:t:B; %массив делений по оси T

z = zeros(21,1); %массив высот

z = min(min(U)):((max(max(U))-min(min(U)))/20):max(max(U));

contour(Y,X,U,z);

colorbar;

title('График изолиний');

grid;

xlabel('Ось T ');

ylabel('Ось X ');

pause;

mesh(U);

colorbar;

title('3D график полученного решения U(x,t)');

xlabel('Ось T ');

ylabel('Ось X ');

zlabel('Ось U');

m=1;

n=1;

for i=1:M/10:M+1

for j=1:N/5:N+1

result(m,n)=U(i,j);

n=n+1;

end

n=1;

m=m+1;

end

result

Текст функции, реализующий решение задачи методом неявной схемы:

clear;

M = input('Введите число узлов по оси X:');

N = input('Введите число узлов по оси T:');

A=1;

B=1/4;

h=A/M; %длина шага разбиения по оси X

t=B/N; %длина шага разбиения по оси Т

g=(t/h)^2;

U=zeros(M+1,N+1); %матрица решений

func=zeros(M+1,N+1);

f=zeros(M+1,1);

a=zeros(M,1);

b=zeros(M,1);

c=zeros(M+1,1);

alpha=zeros(M,1);

beta=zeros(M,1);

x = zeros(M+1,1);

y = zeros(N+1,1);

for i = 1:(M+1),

x(i) = (i-1)*h;

end;

for j = 1:(N+1),

y(j) = (j-1)*t;

end;

p=zeros(M+1,1); %первое начальное условие

q=zeros(M+1,1); %второе начальное условие

fi=zeros(N+1,1); %первое граничное условие

for i=1:M+1

p(i)=100*exp(1-2*x(i));

q(i)=1200*exp(1-2*x(i));

U(i,1)=p(i);

end;

for j=1:N+1

fi(j)=100*exp(12*y(j)+1);

U(1,j)=fi(j);

end;

for i=1:M+1

for j=1:N+1

func(i,j)=14000*exp(12*y(j)-2*x(i)+1);

end;

end;

%нахожнение значений на первом слое

%вычисление коэффициентов для метода прогонки

c(1)=1;

b(1)=0;

a(1)=g/2;

for i=2:M

c(i)=1+g;

a(i)=g/2;

b(i)=g/2;

end;

c(M+1)=1;

a(M)=0;

%вычисление коэффицентов альфа

alpha(1)=b(1)/c(1);

for i=1:M-1

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

end;

%вычисление значений коэф. уравнения f

f(1)=fi(2);

for i=2:M

f(i)=p(i)+t*q(i)+((t^2)/2)*func(i,2);

end;

f(M+1)=g*(p(M)-p(M+1)+h*(-2*U(M+1,1)))+p(M+1)+t*q(M+1)+(t^2/2)*func(M+1,1);

%вычисление коэффицентов бета

beta(1)=f(1)/c(1);

for i=1:M-1

beta(i+1)=(f(i+1)+a(i)*beta(i))/(c(i+1)-a(i)*alpha(i));

end;

U(M+1,2)=(f(M+1)+a(M)*beta(M))/(c(M+1)-a(M)*alpha(M));

%заполнение первого слоя

for i=M:-1:2

U(i,2)=alpha(i)*U(i+1,2)+beta(i);

end;

%нахожнение значений на остальных слоях

%вычисление коэффициентов для метода прогонки

a(1)=g;

for i=2:M

c(i)=1+2*g;

a(i)=g;

b(i)=g;

end;

%вычисление коэффицентов альфа

for i=1:M-1

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

end;

%вычисление значений коэф. уравнения f

for j=3:N+1

f(1)=fi(j);

for i=2:M

f(i)=2*U(i,j-1)-U(i,j-2)+(t^2)*func(i,j-1);

end;

f(M+1)=2*g*(U(M,j-1)-U(M+1,j-1)+h*(-2*U(M+1,j-1)))+func(M+1,j-1)*t^2+2*U(M+1,j-1)-U(M+1,j-2);

%вычисление коэффицентов бета

beta(1)=f(1)/c(1);

for i=1:M-1

beta(i+1)=(f(i+1)+a(i)*beta(i))/(c(i+1)-a(i)*alpha(i));

end;

%вычисление массива решений

U(M+1,j)=f(M+1);

for i=M:-1:2

U(i,j)=alpha(i)*U(i+1,j)+beta(i);

end;

end;

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

y=zeros(M+1,11);

for i=1:M+1

for j=1:11

y(i,j)=U(i,(j-1)*N/10+1);

end;

end;

i=0:1:M;

plot(i/M,y);

legend('t=0','t=1/10 T','t=2/10 T','t=3/10 T','t=4/10 T','t=5/10 T','t=6/10 T','t=7/10 T','t=8/10 T','t=9/10 T','t=T','Location','NorthEastOutside')

title('Параметрический график при различных t');

xlabel('Ось X ');

ylabel('Ось U ');

pause;

X = 0:h:A; %массив делений по оси X

Y = 0:t:B; %массив делений по оси T

z = zeros(21,1); %массив высот

z = min(min(U)):(max(max(U))-min(min(U)))/20:max(max(U));

contour(Y,X,U,z);

colorbar;

title('График изолиний');

grid;

xlabel('Ось T ');

ylabel('Ось X ');

pause;

mesh(U);

colorbar;

title('3D график полученного решения U(x,t)');

xlabel('Ось T ');

ylabel('Ось X ');

zlabel('Ось U');

m=1;

n=1;

for i=1:M/10:M+1

for j=1:N/5:N+1

result(m,n)=U(i,j);

n=n+1;

end

n=1;

m=m+1;

end

result

Текст функции, реализующий решение модельной задачи методом явной схемы:

clear;

t=1;

h=0;

while (t>h)

M = input('Введите число узлов по оси X:');

N = input('Введите число узлов по оси T:');

A=1;

B=1/4;

h=A/M; %длина шага разбиения по оси X

t=B/N; %длина шага разбиения по оси Т

if (t>h) disp('НЕ ВЫПОЛНЕНО УСЛОВИЕ СХОДИМОСТИ'); end;

end;

g=(t/h)^2;

U=zeros(M+1,N+1); %матрица решений

Ur=zeros(M+1,N+1); %модельная функция

f=zeros(M+1,N+1); %матрица значений функции fi,j

x = zeros(M+1,1);

y = zeros(N+1,1);

for i = 1:(M+1),

x(i) = (i-1)*h;

end;

for j = 1:(N+1),

y(j) = (j-1)*t;

end;

p=zeros(M+1,1); %первое начальное условие

q=zeros(M+1,1); %второе начальное условие

fi=zeros(N+1,1); %первое граничное условие

psi=zeros(N+1,1); %второе граничное условие

for i=1:M+1

p(i)=x(i)^5;

q(i)=0;

U(i,1)=p(i);

end;

for j=1:N+1

fi(j)=0;

%psi(j)=y(j)^2*exp(y(j)^2);

U(1,j)=fi(j);

end;

%формирование матрицы значений функции f(i,j)

for i=1:M+1

for j=1:N+1

f(i,j)=-x(i)^3*cos(y(j))*(x(i)^2+20);

Ur(i,j)=x(i)^5*cos(y(j));

end;

end;

%заполнение первого слоя

for i=2:M

U(i,2)=p(i)+t*q(i)+((t^2)/2)*((p(i-1)-2*p(i)+p(i+1))/(h^2)+f(i,1));

end;

U(M+1,2)=g*(p(M)-p(M+1)+h*(-U(M+1,1)+6*cos(y(1)))) +p(M+1)+t*q(M+1)+(t^2/2)*f(M+1,1);

%основной цикл – вычисление массива решений

for j=2:N

for i=2:M

U(i,j+1)=2*(1-g)*U(i,j)+g*(U(i-1,j)+U(i+1,j))-U(i,j-1)+t^2*f(i,j);

end;

U(M+1,j+1)=2*g*(U(M,j)-U(M+1,j)+h*(-U(M+1,j)+6*cos(y(j))))+f(M+1,j)*t^2+2*U(M+1,j)-U(M+1,j-1);

end;

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

y1=zeros(M+1,11);

y2=zeros(M+1,11);

for i=1:M+1

for j=1:11

y1(i,j)=U(i,(j-1)*N/10+1);

y2(i,j)=Ur(i,(j-1)*N/10+1);

end;

end;

i=0:1:M;

subplot(2,1,1)

plot(i/M,y1);

title('Параметрический график при различных t найденной функции');

legend('t=0','t=1/10 T','t=2/10 T','t=3/10 T','t=4/10 T','t=5/10 T','t=6/10 T','t=7/10 T','t=8/10 T','t=9/10 T','t=T','Location','NorthEastOutside')

xlabel('Ось X ');

ylabel('Ось U ');

subplot(2,1,2)

plot(i/M,y2);

legend('t=0','t=1/10 T','t=2/10 T','t=3/10 T','t=4/10 T','t=5/10 T','t=6/10 T','t=7/10 T','t=8/10 T','t=9/10 T','t=T','Location','NorthEastOutside')

title('Параметрический график при различных t модельной функции');

xlabel('Ось X ');

ylabel('Ось U ');

pause;

X = 0:h:A; %массив делений по оси X

Y = 0:t:B; %массив делений по оси T

z = zeros(21,1); %массив высот

z = min(min(U)):((max(max(U))-min(min(U)))/20):max(max(U));

subplot(2,1,1)

contour(Y,X,U,z);

colorbar;

title('График изолиний найденной функции');

grid;

xlabel('Ось T ');

ylabel('Ось X ');

subplot(2,1,2)

contour(Y,X,Ur,z);

colorbar;

title('График изолиний модельной функции');

grid;

xlabel('Ось T ');

ylabel('Ось X ');

pause;

subplot(1,1,1)

mesh(U);

colorbar;

title('3D график найденной функции');

xlabel('Ось T ');

ylabel('Ось X ');

zlabel('Ось U');

pause;

mesh(Ur);

colorbar;

title('3D график модельной функции');

xlabel('Ось T ');

ylabel('Ось X ');

zlabel('Ось U');

maxerror=max(max(abs(Ur-U)));

sprintf('Погрешность вычислений равна %d',maxerror)

Список литературы:

1. Земсков В.Н. Конспект лекций по курсу “Численные методы”.

2. Долголаптев В.Г., Земсков В.Н. Численные методы решения разностных уравнений математической физики. – М., МИЭТ, 1987г.

3. Земсков В.Н., Хахалин С.Я. Метод Сеток. Методические указания к выполнению курсовой работы на персональном компьютере. – М., МИЭТ, 1998.