Конечно-разностный метод решения для уравнений параболического типа (стр. 4 из 4)

аналитическому решению u(x,t) = sin(πx)e^-π2^t+ sin(3πx)e^-9π2^t, истинные значения для последнего представлены в таблице 2

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

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

x_i	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
t₁₁
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), для метода Кранка-Николсона

ЗАКЛЮЧЕНИЕ

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

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

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

Библиографический список

1 Березин И.С., Жидков Н.П. Методы вычислений. Т.2. – М.: Физматгиз, 1962.

2 Мэтьюз, Джон, Г., Финк, Куртис, Д. Численные методы. Использование MATLAB, 3-е издание.— М. : Вильяме, 2001. — 720 с

3 Тихонов А.Н., Самарский А.А. Уравнения математической физики. – М.: Наука, 1972.

4 Формалев В.Ф., Ревизников Д.Л. Численные методы. – М.: ФИЗМАТЛИТ, 2004. - 400 с.

5 Пирумов У.Г. Численные методы. – М.: Издательство МАИ, 1998.

6 Калиткин Н.Н. Численные методы. – М.: Наука, 1976.

ПРИЛОЖЕНИЕ А

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

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

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))]);