Вот собственно код. В результат пишем А и В в координате 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
![$\[\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.\]$](https://dxdy-04.korotkov.co.uk/f/7/b/9/7b93ca317b132c16c8e9b39dbb37b17082.png "$\[\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.\]$")
Одно могу сказать точно. Так как в первом уравнении нелинейная добавка ограничена, то решение