СЛАУ, метод прогонки, проблема накопления ошибки

[pi-314]

Добрый день,
прошу помочь советом по следующей теме:

решаю численно систему УРЧП, в котором одно из уравнений - что-то вроде нелинейного уравнения теплопроводности. При решении СЛАУ для трехточечной по пространству разностной схемы для данного уравнения (нелинейные коэффициенты берутся с предыдущего временного шага) методом прогонки столкнулся с проблемой лавинообразного накопления ошибки при превышении номера шага по времени определенного значения.
Моя первая догадка - значения элементов матрицы очень малы, проблема в машинном округлении чисел. В подтверждение этого, переход от типа Extented к Real - усугубляет проблему и ошибка начинает накапливатьcя еще раньше по времени.
Вторая догадка - источник "тепла" расположен в центре одномерной области задачи, а по краям -ГУ типа постоянное значение - а метод прогонки всегда бегает от одного края к другому, что, наверное, может приводить к численной "асимметрии". В подтверждении этого, ошибка (в виде нефизичных значений искомой переменной) наблюдается только в одной из половинок области задачи.
Буду рад любым комментариям и советам.

[Maik2013]

pi-314
Может систему с граничными условиями и программу подставьте тут для разборка.

[pi-314]

Maik2013 в сообщении #845446 писал(а):

pi-314
Может систему с граничными условиями и программу подставьте тут для разборка.

Ох, я постараюсь освоить LaTeX, но в моей задаче такие длинные выражения!) Боюсь напугать

[pi-314]

Вторую догадку:

pi-314 в сообщении #845432 писал(а):

Вторая догадка - источник "тепла" расположен в центре одномерной области задачи, а по краям -ГУ типа постоянное значение - а метод прогонки всегда бегает от одного края к другому, что, наверное, может приводить к численной "асимметрии". В подтверждении этого, ошибка (в виде нефизичных значений искомой переменной) наблюдается только в одной из половинок области задачи.
Буду рад любым комментариям и советам.

вынужден скорректировать:
Ошибка (в виде нефизичных значений искомой переменной) наблюдается только в одной из половинок области задачи - в первой или второй, что влияет на это - не понятно.

[pi-314]

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

[Maik2013]

pi-314
У Вас примерна такая уравнения был или что то другой

$$ \dfrac{d\theta}{dt}=a_{1}\dfrac{d\,^2\theta}{dx^2}-a_{2}\dfrac{d\theta}{dx}+ t_{\star}k_{0}(1-\theta)\exp\left(-\dfrac{E}{R(T_{0}+\theta(T_{\star}-T_{0}))}\right), $$

[pi-314]

2 Maik2013:
примерно такая.
Вся суть в нелинейных коэффициентах и правой части.
Если распределения значений коэффициентов или функции правой части по ОХ для рассматриваемого временного слоя - несимметричное или с большим размахом значений - решение СЛАУ методом прогонки только в одном направлении OX (например, 0 => 1) будет давать ошибку из-за машинного округления. Для длинных по времени расчетов значение ошибки в итоге может стать сопоставимым характерному значению искомой переменной и уже повлиять на устойчивость самой разностной схемы.

[Maik2013]

pi-314
Уравнения который я добавил тоже хочу решит но не получается и тоже использовал метод прогонки вот программа

Используется синтаксис Delphi

procedure TForm1.Button1Click(Sender: TObject);
var N,M:integer;
Q,k0,E,cp,v0,ro2,T0,c2,Tgor,I,J,k,z,shag,JJ:integer;
L,alfai,betai,beta1,h,tau,kap,a1,a2,lamda1,lamda2,
alfa2,alfa1,ro1,R,eta0,ai,bi,ci,tz:real;


beta, Fi,teta:array of array of Extended;// Двойной массив
alfa:array of Extended;// массив

//Константы
begin
   shag:=0;
   N:=15;
   M:=4000;
   JJ:=10;
   alfa1:=0.5;
   alfa2:=0.5;
   lamda1:=0.084;
   lamda2:=4;
   cp:=1000;
   c2:=660;
   ro1:=0.3;
   v0:=1;
   ro2:=1500;
   Tgor:=1480;
   T0:=315;
   Q:=15130000;
   k0:=500000000;
   E:=126000;
   R:=8.34;
   eta0:=0.077;
   kap:=(alfa1*lamda1+alfa2*lamda2)/(ro1*cp+ro2*c2);
   L:=0.1;
   tz:=1;//sqr(L)/kap;
   h:=L/N;
   tau:=0.1;
   a1:= (tz*kap)/sqr(L);
   a2:= (tz*v0)/L;
   ai:=a2/h-a1/sqr(h);
   ci:=-a1/sqr(h);
   bi:=(-2*a1)/(sqr(h))-1/tau+a2/h;
   //SH:=10;

    Memo2.Lines.Add('KAP='+FloatToStr(KAP));
    Memo2.Lines.Add('L='+FloatToStr(L));
    Memo2.Lines.Add('H='+FloatToStr(H));
    Memo2.Lines.Add('A1='+FloatToStr(A1));
    Memo2.Lines.Add('A2='+FloatToStr(A2));
    Memo2.Lines.Add('Ai='+FloatToStr(AI));
    Memo2.Lines.Add('Ci='+FloatToStr(CI));
    Memo2.Lines.Add('Bi='+FloatToStr(BI));
    Memo2.Lines.Add('h='+FloatToStr(h));
    Memo2.Lines.Add('tau='+FloatToStr(tau));
 SetLength(TETA, N+1,M+1);   //Дарозии массыв
 SetLength(Fi, N+1,M+1);
 SetLength(alfa, N+1);
 SetLength(beta, N+1,M+1);
 alfa[1]:=0;
 beta[1,0]:=1;
   for k:=0 to M do
   begin
 teta[0,k]:=1;
   end;
  for z:=1 to N-1 do
  Begin
   teta[z,0]:=0;
  end;
  for J:=0 to M-1 do
begin
      for i:=1 to N-1 do
      begin
   Fi[i,j]:=(-teta[i,j]/tau)-tz*k0*(1-teta[i,j])*exp(-(E)/(R*(T0+teta[i,j]*(Tgor-T0))));
      end;
  for i:=2 to N do
   begin
   alfa[i]:=ai/(bi-ci*alfa[i-1]);
   beta[i,j]:=-(ci*beta[i-1,j]-Fi[i-1,j])/(bi-ci*alfa[i-1]);
   end;
   teta[N,j+1]:=beta[N,j]/(1-alfa[N]);
   for i:=N-1 downto 1 do
   begin
    teta[i,j+1]:=alfa[i+1]*teta[i+1,j+1]+beta[i+1,j];
   end;
   end;
   J:=1;
   repeat
     for i:=0 to N-1 do
    begin
     shag:=shag+1;
     if shag>=1 then
     begin
       //Memo1.Lines.Add(FloatToStr(teta[N,j+1]));
      memo1.Lines.Add('TETA('+INTTOSTR(I)+','+INTTOSTR(J)+') = '+FLOATTOSTR(TETA[I,J]));
       shag:=0;
     end;
    end;
    J:=J+JJ;
 until j>M-1;
end;
end;
 

Я должен получит примерна такой график

но не как не получается. Может Вы какой ни буд советь даете ?