Model_graph.m

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

%Параметры N1, N2 - количество узлов по осям x и t

function [U,U_n] = fun(N1,N2)

%Задание границ области дискретного изменения аргумента

L_const=1;

T_const=2^(3/2)/3;

x_begin=0;

x_end=L_const;

t_begin=0;

t_end=T_const;

%Определение шага по осям x и t

h=L_const/N1;

tau=T_const/N2;

gamma=tau/h;

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

X=x_begin:h:x_end;

T=t_begin:tau:t_end;

N=length(X);

M=length(T);

%Сеточная функция внешней силы f(x,t)

for i=1:N

for j=1:M

f(i,j)=-81/4*X(i)-27/4;

end;

end;

%------------------------------------------------------------------------------------------

%--------------Сеточный аналог теоретического решения----------------------

%------------------------------------------------------------------------------------------

for i=1:N

for j=1:M

U_a(i,j)=27/8*X(i)^3-27/8*T(j)^2;

end;

end;

figure;

mesh(T,X,U_a);

title('Teoretic function');

xlabel('t');

ylabel('x');

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

figure;

contour(U_a,20);

%---------------------------------------------------------------------------------------------------

%-------Решение краевой задачи с применением явной разностной схемы--------

%---------------------------------------------------------------------------------------------------

%Сеточная функция U(x,t)

U=zeros(N,M);

%Задание первого начального условия U(x,0)

for i=1:N

U(i,1)=27/8*X(i)^3;

end;

%Задание первого граничного условия U(L,t)

for j=1:M

U(N,j)=27/8-27/8*T(j)^2;

end

%Задание второго начального условия Ut(x,0) - по известным U(i,0)

%определяем значения U(i,1) для i=1,...,N-1

for i=2:N-1

U(i,2)=(1-gamma^2)*U(i,1)+gamma^2*U(i-1,1)/2+gamma^2*U(i+1,1)/2+tau^2*f(i,1)/2;

end

%Определение значения U(0,1) из начальных и граничных условий

U(1,2)=gamma^2*U(2,1)-(gamma^2-1)*U(1,1)+tau^2*f(1,1)/2;

%Вычисление значения функции в узлах

for j=2:M-1

%Из второго граничного условия находим U(0,j) для j=2,...,M

U(1,j+1)=2*gamma^2*U(2,j)+2*(1-gamma^2)*U(1,j)-U(1,j-1)+tau^2*f(1,j);

%Находим значения сеточной функции во внутренних узлах сетки

for i=2:N-1

U(i,j+1)=2*(1-gamma^2)*U(i,j)+gamma^2*(U(i-1,j)+U(i+1,j))-U(i,j-1)+tau^2*f(i,j);

end

end

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

figure;

mesh(T,X,U);

title('Evident scheme');

xlabel('t');

ylabel('x');

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

figure;

contour(U,20);

%---------------------------------------------------------------------------------------------------

%------Решение краевой задачи с применением неявной разностной схемы-------

%----------------------------------------------------------------------------------------------------

%Сеточная функция U(x,t)

U_n=zeros(N,M);

%Задание первого начального условия U(x,0)

for i=1:N

U_n(i,1)=27/8*X(i)^3;

end;

%Задание первого граничного условия U(L,t)

for j=1:M

U_n(N,j)=27/8-27/8*T(j)^2;

end

%Вычисление значений сеточной функции U(i,1) методом прогонки

%Определение прогоночных коэффициентов

A=zeros(N-1,1);

B=zeros(N-1,1);

ALPHA=zeros(N-1,1);

F=zeros(N,1);

C=zeros(N,1);

BETA=zeros(N,1);

for i=1:N-2

A(i)=gamma^2/2;

B(i+1)=gamma^2/2;

C(i+1)=1+gamma^2;

F(i+1)=U_n(i+1,1)+tau^2*f(i+1,2)/2;

end;

A(N-1)=0;

B(1)=gamma^2;

C(1)=1+gamma^2;

C(N)=1;

F(1)=U_n(1,1)+tau^2*f(1,1)/2;

F(N)=27/8-27/8*T(1)^2;

%Прямой ход прогонки

ALPHA(1)=B(1)/C(1);

BETA(1)=F(1)/C(1);

for i=2:N-1

ALPHA(i)=B(i)/(C(i)-A(i-1)*ALPHA(i-1));

BETA(i)=(F(i)+A(i-1)*BETA(i-1))/(C(i)-A(i-1)*ALPHA(i-1));

end;

BETA(N)=(F(N)+A(N-1)*BETA(N-1))/(C(N)-A(N-1)*ALPHA(N-1));

%Обратный ход прогонки

U_n(N,2)=BETA(N);

for i=N-1:-1:1

U_n(i,2)=ALPHA(i)*U_n(i+1,2)+BETA(i);

end;

%Вычисление значений сеточной функции U(i,j) j=2,...,M методом прогонки

%Определение прогоночных коэффициентов

A=zeros(N-1,1);

B=zeros(N-1,1);

ALPHA=zeros(N-1,1);

F=zeros(N,1);

C=zeros(N,1);

BETA=zeros(N,1);

for i=1:N-2

A(i)=gamma^2;

B(i+1)=gamma^2;

C(i+1)=1+2*gamma^2;

end;

A(N-1)=0;

B(1)=2*gamma^2;

C(1)=1+2*gamma^2;

C(N)=1;

for j=3:M

F=zeros(N,1);

F(1)=2*U(1,j-1)-U(1,j-2)+tau^2*f(1,j-1);

F(N)=27/8-27/8*T(j)^2;

ALPHA=zeros(N-1,1);

BETA=zeros(N,1);

%Прямой ход прогонки

ALPHA(1)=B(1)/C(1);

BETA(1)=F(1)/C(1);

for i=2:N-1

F(i)=2*U(i,j-1)-U(i,j-2)+tau^2*f(i,j-1);

ALPHA(i)=B(i)/(C(i)-A(i-1)*ALPHA(i-1));

BETA(i)=(F(i)+A(i-1)*BETA(i-1))/(C(i)-A(i-1)*ALPHA(i-1));

end;

BETA(N)=(F(N)+A(N-1)*BETA(N-1))/(C(N)-A(N-1)*ALPHA(N-1));

%Обратный ход прогонки

U_n(N,j)=BETA(N);

for i=N-1:-1:1

U_n(i,j)=ALPHA(i)*U_n(i+1,j)+BETA(i);

end;

end;

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

figure;

mesh(T,X,U_n);

title('Not evident scheme');

xlabel('t');

ylabel('x');

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

figure;

contour(U_n,20);