%Функция расчета сеточной функции с использованием явной и неявной схем
%Параметры 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)=-9*(27/4*X(i)^2-3/2^0.5*T(j)-3)/((9/4*X(i)^2+3/2^0.5*T(j)+1)^3);
end;
end;
%---------------------------------------------------------------------------------------------------
%-------Решение краевой задачи с применением явной разностной схемы--------
%---------------------------------------------------------------------------------------------------
%Сеточная функция U(x,t)
U=zeros(N,M);
%Задание первого начального условия U(x,0)
for i=1:N
U(i,1)=2/((3*X(i)/2)^2+1);
end;
%Задание первого граничного условия U(L,t)
for j=1:M
U(N,j)=8/(13+6*2^0.5*T(j));
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)-4*tau/((9*(X(i)^2)/4+1)^2)+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)-4*tau+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)=2/((3*X(i)/2)^2+1);
end;
%Задание первого граничного условия U(L,t)
for j=1:M
U_n(N,j)=8/(13+6*2^0.5*T(j));
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)-4*tau/((9*(((i)*h)^2)/4+1)^2)+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-4*tau;
F(N)=8/(13+(6*2^0.5)*tau);
%Прямой ход прогонки
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)=8/(13+(6*2^0.5)*tau*(j-1));
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);
disp('Построение двумерных параметрических графиков');
%Построение двумерного параметрического графика U(x) при различных t
figure;
hold;
title('U(x), t=const');
xlabel('x');
ylabel('U(x)');
for j=1:M
if mod(j,5)==0
plot(X,U(:,j),'b');
plot(X,U_n(:,j),'r');
end;
end
legend('evident scheme','not evident scheme');
%Построение двумерного параметрического графика U(t) при различных x
figure;
hold;
title('U(t), x=const');
xlabel('t');
ylabel('U(t)');
for i=1:N
if mod(i,5)==0
plot(X,U(i,:),'b');
plot(X,U_n(i,:),'r');
end;
end
legend('evident scheme','not evident scheme');