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

 

 

 

 

 

 

 

 

 

 

 

3. Расчетная часть.

3.1. Математическое описание решения задачи

С учетом выше изложенной теории:

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

 

 

 

 

 

 

 

 

2) для схемы Кранка-Николсона преобразовываем уравнение

 к  виду:

 

 

 

где

 

 

 

3.2. Алгоритм решения задачи

 

 

На основе выше изложенной теории, составляем блок-схему явной разностной схемы.

 

Описание алгоритма явной разностной системы.

 

1. Начало алгоритма.
2. Задать начальные краевые условия задачи: _ф (_Х) := (t) := 0 (t) := О

 

Перейти к следующему шагу.

3. Ввести N - число точек сетки по х, К - число точек сетки по t, L - длина отрезка по х, Т - длина отрезка по t. Перейти к следующему шагу.
4. Производится расчет шагов по х и t соответственно: h = L/N, ∆ = Т/К. Перейти к следующему шагу.
5. Цикл. Изменять i от 0 до N. Для i (от 0 до N)  перейти к следующему шагу, 
   для i>N перейти к шагу 8. , . ,
6. Присвоить значение узла сетки Х[ := i*h Перейти-к следующему шагу.
7. Вычисляются массивы начальных условий: ui:=:.ф(Xi)
8. Цикл. Изменять j от 0 до К. Для j (от 0 до К) перейти к следующему шагу, для j>К перейти к шагу 10.
9. Присвоить значение узла сетки tj:=j*∆. Перейти к следующему шагу.

 

10. Вычисляется коэффициент у=∆ /h². Перейти к следующему шагу.
11. Цикл. Изменять j от 0 до К. Для j (от 0 до К) перейти к следующему шагу 14.
12. Вычисляются массивы условий: U0j:= µ(tj). Перейти к следующему шагу.
13. Вычисляются массивы краевых условий: Unj:=v(tj). Перейти к следующему шагу.
14. Цикл. Изменять j от 0 до К-1. Для j (от 0 до К-1) перейти к следующему шагу, для j>К-1 перейти к шагу 17.
15. Цикл. Изменять i от 0 до N-1. Для i (от 0 до N-1) перейти к следующему шагу, для i>N-1 перейти к шагу 14.
16. Поиск окончательного решения в узлах сетки: Uij:= у*Ui-1j + (1 -2* у)*Uij + y*Ui+1j.Перейти к следующему шагу.
17. Цикл. Изменять i от 0 до N. Для i (от 0 до N) перейти к следующему шагу, для i >N перейти к шагу 20.
18. Цикл. Изменять j от 0 до К. Для j (от 0 до К) перейти к следующему шагу, для j>К перейти к шагу 17.
19. Вычисляются точные значения функции в узлах сетки и ukij=E(xi,tj), Rij:= i(Uij-Ukij)/Ukij 100i. Оценивается относительная погрешность как модуль отношения разности между точными значениями и расчетными на точные значения.
20. Вывести вектор u,R.

 

3.3. Решение задачи в MatLab.

 

   Решение будем искать в ППП MatLab 7. Создадим четыре выполняемых m-фала: crnich.m – файл-функция с реализацией метода Кранка-Николсона; trisys.m – файл-функция метода прогонки; f.m – файл-функция задающая начальное условие задачи; fе.m – файл-функция задающая функцию определяющую точное решение задачи(найдена аналитическим путем).

Для простоты возьмем шаг Δх = h = 0,1. Таким образом, соотношение

 L =1. Пусть решетка имеет n=10 столбцов в ширину.

Анализ  результатов

Решения для данных параметров отразим в  таблице 1. Трехмерное изображение данных из таблицы покажем на рисунке 5.

 

Таблица 1 – Значения u(х_i, t_i), полученные методом Кранка-Николсона

xi	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
ti
0	0	1.1180	1.5388	1.1180	0.3633	0	0.3633	1.1180	1.5388	1.1180	0
0.01	0	0.6169	0.9288	0.8621	0.6177	0.4905	0.6177	0.8621	0.9288	0.6169	0
0.02	0	0.3942	0.6480	0.7186	0.6800	0.6488	0.6800	0.7186	0.6480	0.3942	0
0.03	0	0.2887	0.5067	0.6253	0.6665	0.6733	0.6665	0.6253	0.5067	0.2887	0
0.04	0	0.2331	0.4258	0.5560	0.6251	0.6458	0.6251	0.5560	0.4258	0.2331	0
0.05	0	0.1995	0.3720	0.4996	0.5754	0.6002	0.5754	0.4996	0.3720	0.1995	0
0.06	0	0.1759	0.3315	0.4511	0.5253	0.5504	0.5253	0.4511	0.3315	0.1759	0
0.07	0	0.1574	0.2981	0.4082	0.4778	0.5015	0.4778	0.4082	0.2981	0.1574	0
0.08	0	0.1419	0.2693	0.3698	0.4338	0.4558	0.4338	0.3698	0.2697	0.1419	0
0.09	0	0.183	0.2437	0.3351	0.3936	0.4137	0.3936	0.3351	0.2437	0.1283	0
0.1	0	0.1161	0.2208	0.3038	0.3570	0.3753	0.3570	0.3038	0.2208	0.1161	0

 

Величины, полученные методом Кранка-Николсона, достаточно близки к

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

Максимальная  погрешность для данных параметров равна 0,005

 

Таблица 2 – точные значения u(х_i, t_i), при t=0.1

xi	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
t11
0.1	0	0.1153	0.2192	0.3016	0.3544	0.3726	0.3544	0.3016	0.2192	0.1153	0

 

Рисунок 5 – Решение u=u(х_i, t_i), для метода Кранка-Николсона

 

 

Программа для расчета по методу Кранка-Николсона

 

f.m

function F=f(x)

F=sin(pi*x)+sin(3*pi*x);

 

ft.m

function F=f(x,t)

F=sin(pi*x)*exp(-pi^2*t)+sin(3*pi*x)*exp(-9*pi^2*t);

 

tisys.m

function X=trisys(A,D,C,B)

N=length(B);

for k=2:N

mult=A(k-1)/D(k-1);

D(k)=D(k)-mult*C(k-1);

B(k)=B(k)-mult*B(k-1);

end

X(N)=B(N)/D(N);

for k= N-1:-1:1

X(k)=(B(k)-C(k)*X(k+1))/D(k);

end

 

crnich.m

 

function [U,Y]=crnich(c1,c2,a,b,c,n,m)

clc;

% - c1=u(0,t) и c2=u(a,t)

% - а и b - правые точки интервалов [0,а] и [0,Ь] 

% - с - постоянная уравнения теплопроводности

% - n и m - число точек решетки  на интервалах [0,а] и [0,Ь] 

%Выход  - U - матрица решений 

 

%Инициализация  параметров и матрицы U

h=a/(n-1);

k=b/(m-1);

r=c^2*k/h^2;

s1=2+2/r;

s2=2/r-2;

U=zeros(n,m);

 

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

U(1,1:m)=c1;

U(n,1:m)=c2;

 

%Генерирование  первого ряда 

U(2:n-1,1)=f(h:h:(n-2)*h)';

 

%Формирование  диагональных и не лежащих на диагонали

%элементов А, вектора постоянных В '

%и  решение трехдиагональной системы АХ=В

Vd(1,1:n)=s1*ones(1,n);

Vd(1)=1;

Vd(n)=1;

Va=-ones(1,n-1);

Va(n-1)=0;

Vc=-ones(1,n-1);

Vc(1)=0;

Vb(1)=c1;

Vb(n)=c2;

for j=2:m

for i=2:n-1

Vb(i)=U(i-1,j-1)+U(i+1,j-1)+s2*U(i,j-1);

end

X=trisys(Va,Vd,Vc,Vb);

U(1:n,j)=X';

end

U=U';

%точное  решение и определение погрешности

x=0:h:h*(n-1);

t=0:k:k*(m-1);

for i=1:1:n

for j=1:1:m

Y(i,j)=ft(x(i),t(j));

end

end

Y=Y';

pogr=(abs(Y-U));

max(max(pogr))

surf(U)

xlabel('X');

ylabel('t');

zlabel('U(x,t)');

colorbar;

axis([0 n 0 m 0 max(max(U))]);

 

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

 

ЗАКЛЮЧЕНИЕ

 

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

В МКР строятся, как правило, регулярные сетки, особенности геометрии области  учитываются только в около граничных  узлах. В связи с этим МКР чаще применяется для анализа задач  с прямолинейными границами областей определения функций.

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

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Список используемой литературы:

1. Амосов А.А, Дубянский Ю.А., Копченов Н.В. Вычислительные методы для инженеров: Учеб. пособие. –  М.: Высш. шк., 2005.–387
2. Дьяконов В.П. MathCAD 2001. специальный справочник. – СПб.: Питер, 2004. –583
3. Заварыкин В.М. Численные методы: учеб. пособие для студентов пед. ин-в. – М.: Просвещение, 2006. – 485
4. Кирьянов Д.В. Самоучитель MathCAD 2001. – СПб.: БХВ – Петербург, 2007. – 539
5. Макаров Е.Г. Инженерные расчеты в MathCAD. Учебный курс. – СПб.: Питер, 2009. – 528.