Численное решение уравнений теплопроводности методом разностных схем

%l - длинна стержня

%N - Разбиение сетки по пространственной координате

%T - Разбиение сетки по временной координате

%k_U - функция k(U)

%f_U - функция f(U)

%U0 - Начальное распределение температуры

%U1 - Функция подачи тепла на начало стержня

%U2 - Функция подачи тепла на конец стержня

%E - точность вычислений

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

%создание сетки

h=l/N;

tau=T/K;

x=[0:h:h*N];

t=[0:tau:tau*K];

coef=0.5*tau/h^2;

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

Ymain=zeros(K+1,N+1); % матрица полученных значений.

Ymain(1,:)=feval(U0,x); % Вставка значений начальных условий в первую строку, где U(x,0)=f(x);

Ymain(:,1)=feval(U1,t'); % Вставка значений начальных условий в первый столбец, где U(0,t)=f1(t);

Ymain(:,N+1)=feval(U2,t'); % Вставка значений начальных условий в последний столбец, где U(1,t)=f2(t);

U=zeros(N-1,N-1);

D=1;

%0-я итерация  для заполнения значений для  1-й итерации...

for i=1:K

    %заполняем  ai

    for j=1:N

         a(j)=0.5*(feval(k_U,Ymain(i,j+1))+feval(k_U,Ymain(i,j))); %Всего  N элементов, крайние подут в  правые части

    end;

    %составление  основной матрицы для С.Л.А.У.

    %диагонали, расположенные вышеи ниже главной

    UpDiag=(tau/h^2).*a;

    U=diag(UpDiag(1,2:N-1),1)+diag(UpDiag(1,2:N-1),-1);

    %Главная диагональ

    for j=1:N-1

        U(j,j)=-(1+(tau/h^2)*(a(j)+a(j+1)));

    end;

    % Составление  правых частей для С.Л.А.У.

    B=-tau*feval(f_U,Ymain(i,2:N))-Ymain(i,2:N);

    B(1,1)=B(1,1)-(tau/h^2)*a(1)*Ymain(i+1,1);

    B(1,N-1)=B(1,N-1)-(tau/h^2)*a(N)*Ymain(i+1,N+1);

    Solution=U\B';

    Ymain(i+1,2:N)=Solution';

end;

%Далее запускаем  итерационный процесс

while(E<D)

    for i=1:K

       %заполняем ai

       for j=1:N

         a(j)=0.5*(feval(k_U,Ymain(i+1,j+1))+feval(k_U,Ymain(i+1,j))); %Всего  N элементов, крайние подут в  правые части

       end;

       %составление основной матрицы  для С.Л.А.У.

       %диагонали, расположенные вышеи  ниже главной

       UpDiag=(tau/h^2).*a;

       U=diag(UpDiag(1,2:N-1),1)+diag(UpDiag(1,2:N-1),-1);

       %Главная диагональ

       for j=1:N-1

          U(j,j)=-(1+(tau/h^2)*(a(j)+a(j+1)));

       end;

       % Составление правых частей для  С.Л.А.У.

       B=-tau*feval(f_U,Ymain(i+1,2:N))-Ymain(i,2:N);

       B(1,1)=B(1,1)-(tau/h^2)*a(1)*Ymain(i+1,1);

       B(1,N-1)=B(1,N-1)-(tau/h^2)*a(N)*Ymain(i+1,N+1);

       Solution=U\B';

       M=norm((Ymain(i+1,2:N)-Solution'));

       %D=norm((Ymain(i+1,2:N)-Solution'));

       if(M<D)

          D=M;

       end;

       Ymain(i+1,2:N)=Solution';

    end;

end;

figure;

surf(x,t,Ymain);

   %ansK_U=Y^2;

   %ansK_U=1;

   %ansK_U=Y.^7; %степенная  нелинейность

   ansK_U=-exp(3*Y); %экспоненциальная нелинейность

   %Y=0; %без источника

   %Y=sin(U);

   Y=0.45*U.*(1.-U); %Уравнение Фишера.

   %Y=0.45*U-0.1*U.^3; %Уравнение Ньюэлла--Уайтхеда (Newell--Whitehead).

   UX=sin(X+0.45);

  UT=cos(T+0.45);

  UT=6*T+0.9;