Решение нелинейной системы дифференциальных уравнений.

[tereshin]

Здравствуйте!
Есть следующая система уравнений.

$\[\left\{ \begin{array}{l} \frac{{\partial A}}{{\partial t}} = ( - \sin 2\pi z)\frac{B}{{1 + {B^2} + {{\left( {\frac{{\partial A}}{{\partial z}}} \right)}^2}}} + \frac{{{\partial ^2}A}}{{\partial {z^2}}}\\ \frac{{\partial B}}{{\partial t}} = \frac{\partial }{{\partial z}}\left( { - \sin 2\pi z\frac{{\frac{{\partial A}}{{\partial z}}}}{{1 + {B^2} + {{\left( {\frac{{\partial A}}{{\partial z}}} \right)}^2}}}} \right) - \frac{{\partial A}}{{\partial z}} + \frac{{{\partial ^2}B}}{{\partial {z^2}}} \end{array} \right.\]$

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

[tereshin]

Вот собственно код. В результат пишем А и В в координате z=0,3 при времени от 0 до 50.

program Potsdam_dynamo

implicit none 
real*8, allocatable ::  A_1(:),B_1(:),A_2(:),B_2(:),z(:)
real*8 pi/3.141592653589793/

integer n,ig,i,nom
real*8 C_omega,C_alpha,dz,dt,t
real*8 centr_diff,right_diff1,right_diff2,diff_2nd_A,diff_2nd_B
real*8 znam0,drob1,drob2,drob3,drob4
!---------------------Task Parameters section--------------------------------
n=101 !number of nodes
C_omega = 50.0d0
C_alpha = 20.0d0

dz = 1.0d0/(n-1)
dt=dz/300.0d0   !convergence condition!!!!!
!----------------------Memory Allocation for Arrays--------------------------
allocate(A_1(n))
allocate(A_2(n))
allocate(B_1(n))
allocate(B_2(n))
allocate(z(n)) 
!Open Files for future results writing
open(1001,file="resultA.txt")
open(1002,file="resultB.txt")
!----------------------Initializing Grid and Time----------------------------
do i=1,n
  z(i) = (i-1) * dz
enddo
t=0.0d0
ig=0
!----------------------Initial Condition Section-----------------------------
B_1=sin(2*pi*z)
A_1(1:n)=0.0d0

B_2(1:n) = 0.0d0
A_2(1:n) = 0.0d0
!--------------------Main Loop (Time loop)------------------------------------
do while (t<50.0d0)
    ig=ig+1
      t=ig*dt
    !-------------------Spatio Loop-------------------------------------------
      do i=2,n-1
        !Differences
            centr_diff = (A_1(i+1) - A_1(i-1))/(2*dz);
        right_diff1 = (A_1(i+1) - A_1(i))/dz;
        right_diff2 = (A_1(i) - A_1(i-1))/dz;
        diff_2nd_A = (A_1(i+1) - 2.0d0*A_1(i) + A_1(i-1))/(dz*dz);
        diff_2nd_B = (B_1(i+1) - 2.0d0*B_1(i) + B_1(i-1))/(dz*dz);
        !for convenience
            znam0 = 1.0d0 + B_1(i)*B_1(i);
        drob1 = B_1(i)/(znam0+centr_diff*centr_diff);
        drob2 = centr_diff/(znam0+centr_diff*centr_diff);
        drob3 = right_diff1/(znam0+right_diff1*right_diff1);
        drob4 = right_diff2/(znam0+right_diff2*right_diff2);
        !Finite Difference Scheme Realization
        A_2(i) = (C_alpha*(-sin(2*pi*z(i)))*drob1+diff_2nd_A)*dt + A_1(i);
        B_2(i) = (-C_alpha* ( -2*pi*cos(2*pi*z(i)) *drob2 -sin(2*pi*z(i)) * ((drob3-drob4)/dz) ) - C_omega*centr_diff + diff_2nd_B)*dt+B_1(i);
    enddo
      !------------------------Boundary Condition Section-------------------------
    A_2(1) = A_2(2);
    A_2(n) = A_2(n-1);
    B_2(1) = 0.0d0;
    B_2(n) = 0.0d0;
    !-------------Preparation for the transition to the next time step----------    
    A_1(1:n) = A_2(1:n);
    B_1(1:n) = B_2(1:n);
      !-------------------------Printing Results to File--------------------------
    if (mod(ig,40)==0) then
   nom = 0.3*n;
        write(1001,*)A_1(nom);
        write(1002,*)B_1(nom);
   print *,'Time = ',t
   endif
    
enddo

!-------------------Deallocating Memory------------------------------------------
deallocate(A_1)
deallocate(A_2)
deallocate(B_1)
deallocate(B_2)
deallocate(z) 
!-------------------Closing Files------------------------------------------------
close(1001)
close(1002)

end program Potsdam_dynamo

[sup]

С теоретической точки зрения, система вполне приличная. Поэтому скорее всего, у Вас ошибки в расчетах. Правда, Вы забыли добавить начальные и краевые условия. Судя по коду, на границе функции равны 0, а при $t=0$ начальные данные $A \equiv 0, B=\sin 2\pi z$.
Я быстренько попробовал, и у меня вроде все "устаканивается". В Вашем коде присутствуют некие "магические" константы C_alpha, C_omega.
Откуда они взялись? Вот, например,

A_2(i) = (C_alpha*(-sin(2*pi*z(i)))*drob1+diff_2nd_A)*dt + A_1(i);


В этой строчке вроде бы никакой C_alpha быть не должно.

[tereshin]

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

[sup]

Константы и вправду на сходимость не влияют. Но все же странно, зачем нужна улучшенная аппроксимация на краю. Может проблема просто "спряталась"?
Как я понял, вы во втором уравнении честно продифференцировали и все слагаемые надлежащим образом записали. А почему бы Вам не использовать разность вместо дифференцирования? Я, лично, так и сделал. Более того, я бы все таки применил неявную схему для линейных слагаемых.

[tereshin]

Улучшенная аппроксимация на краю была нужна для того чтобы схема получила второй порядок аппроксимации по пространству и начала быстрее сходится. Оказывается она и раньше сходилась у меня но я не мог это наблюдать ибо плохо учтенные краевые условия давали огромную ошибку с течением времени. В качестве результата я должен был выдавать изменение решения в одной точке с течением времени (время до 50) и как вы понимаете там сильно все портилось.
Спасибо огромное за последний совет насчет использования неявной схемы с линейными членами. Постараюсь воплотить его в жизнь.

[sup]

Дело в том, что "нелинейщина" дает не такой уж и "большой" вклад, поэтому её можно брать с предыдущего слоя. А вот неявная схема гарантирует наличие хороших оценок на разностное решение. Это в принципе позволит Вам уменьшить разбиение по $t$. Можно даже и так. Сначала решить первое уравнение неявной схемой. А потом во второе уже подставлять найденное $A$. Но это так ... мысли вслух.

[tereshin]

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

[sup]

А я это так ботнул, без особой ответственности. Одно могу сказать точно. Так как в первом уравнении нелинейная добавка ограничена, то решение $A$ обладает массой "хороших" свойств. Во втором уравнении нелинейная добавка - это пространственная производная от ограниченной функции. Тоже неплохо. Все это позволяет "хулиганить" в свое удовольствие. Думаю, результат от этого не слишком изменится. Можно попробовать и так и так ...и использовать тот подход, который быстрее сходится. В общем случае такие шалости, скорее всего, не пройдут. Но данная система уж очень специфическая.