Температурный расчет с помощью вычислений информационной математики
c-----------------------------------------------------------------
cПОДПРОГРАММА РЕШЕНИЯ СИСТЕМЫ УРАВНЕНИЙ МЕТОДОМ ЗЕЙДЕЛЯ
c
с integer N-входное количество уравнений
c real y(6,N)-входной массив уравнений,содержащий следующие поля:
c y(1,N)-номер точки по оси X
c y(2,N)-номер точки по оси Y
c y(3,N)-коэфициен уравнения для Q(y(1,N)-1,y(2,N))
c y(3,N)=h^2/(2*(h^2+k^2))
c y(4,N)-коэфициен уравнения для Q(y(1,N),y(2,N)-1)
c y(4,N)=k^2/(2*(h^2+k^2))
c y(5,N)-коэфициен уравнения для Q(y(1,N)+1,y(2,N))
c y(5,N)=h^2/(2*(h^2+k^2))
c y(6,N)-коэфициен уравнения для Q(y(1,N),y(2,N)+1)
c y(6,N)=k^2/(2*(h^2+k^2))
c integer M-число узлов по оси X
c integer P-число узлов по оси Y
c real Q(M,P)-входной массив начальных значений Y
c real Q(M,P)-выходной массив вычисленых значений Y
c real E-погрешность вычислений
c------------------------------------------------------------------
subroutine zeidel(N,y,M,P,q,E)
integer N,M,P,I,S
real y(6,N),q(M,P),E,EI,NEXTQ
c------------------------------------------------------------------
c вычисление коэфициента сходимости процесса
c mj=y(5,1)+y(6,1)
c и выражения
c km=mj/(1-mj)
C НО Т.К. MJ=0.5 ТО KM=1 И СЛЕДОВАТЕЛЬНО ЕГО МОЖНО ОПУСТИТЬ
c-----------------------------------------------------------------
c KM=(y(5,1)+y(6,1))/(1-y(5,1)+y(6,1))
c------------------------------------------------------------------
c итерации
c S-счетчик итераций
c------------------------------------------------------------------
S=0
1 EI=0.
S=S+1
do I=1,N
NEXTQ=y(3,i)*Q(y(1,i)-1,y(2,i))+
+ y(4,i)*Q(y(1,i),y(2,i)-1)+
+ y(5,i)*Q(y(1,i)+1,y(2,i))+
+ y(6,i)*Q(y(1,i),y(2,i)+1)
c------------------------------------------------------------------
c вычисление погрешности на данной итерации
c------------------------------------------------------------------
if (abs(NEXTQ-q(y(1,i),y(2,i))).gt.EI)
+ EI=abs(NEXTQ-q(y(1,i),y(2,i)))
c print *,'x=',y(1,i),' y=',y(2,i)
q(y(1,i),y(2,i))=NEXTQ
enddo
c print '(16h Итерация номер ,i5,13h погрешность=,E15.7)',S,EI
if (EI.gt.E)goto 1
c do i=P,1,-1
c print '(21e10.3)',(q(j,i),j=1,M)
c enddo
end
ТЕСТ
В качестве теста выполним одну итерацию для системы , полученной в предыдущем пункте.
при начальных условиях:
все остальные Yij:=0.
Получается:
Результат:
Полный текст программы.
c------------------------------------------------------------------
c ПОДПРОГРАММА СОСТАВЛЕНИЯ СИСТЕМЫ УРАВНЕНИЙ
c МЕТОДОМ КОНЕЧНЫХ РАЗНОСТЕЙ
c
c real H-шаг по оси X
c real K-шаг по оси Y
c real N-количество уравнений(примерное число,желательно N=M*P)
c real y(6,N)-выходной массив уравнений,содержащий следующие поля:
c y(1,N)-номер точки по оси X
c y(2,N)-номер точки по оси Y
c y(3,N)-коэфициен уравнения для Q(y(1,N)-1,y(2,N))
c y(3,N)=h^2/(2*(h^2+k^2))
c y(4,N)-коэфициен уравнения для Q(y(1,N),y(2,N)-1)
c y(4,N)=k^2/(2*(h^2+k^2))
c y(5,N)-коэфициен уравнения для Q(y(1,N)+1,y(2,N))
c y(5,N)=h^2/(2*(h^2+k^2))
c y(6,N)-коэфициен уравнения для Q(y(1,N),y(2,N)+1)
c y(6,N)=k^2/(2*(h^2+k^2))
c integer M-число узлов по оси X
c integer P-число узлов по оси Y
c real Q(M,P)-массив значений Y
c integer N-выходное количество получившихся уравнений
c------------------------------------------------------------------
subroutine mkr(H,K,N,y,M,P,q)
integer M,P,IIX,IIY,NN,N,KR1,KR2,KR3
real y(6,N),H,K,q(M,P),HX,KY
c-----------------------------------------------------------------
c подсчитываю коэфициенты
c h^2/(2*(h^2+k^2))
c и
c k^2/(2*(h^2+k^2))
c-----------------------------------------------------------------
HX=H**2/(2*(H**2+K**2))
KY=K**2/(2*(H**2+K**2))
c-----------------------------------------------------------------
c составление уравнений
c и
c присваивание начальных значений
c
c nn-счетчик уровнений
c iix-номер текущего узла по оси X
c iiy-номер текущего узла по оси Y
c-----------------------------------------------------------------
NN=0
KR1=((P-1)/8)*3+1
KR2=((P-1)/8)*5+1
KR3=((M-1)/4)*3+1
do IIY=2,P-1
do IIX=2,M
if (NN.eq.N)then
print *,'ПЕРЕПОЛНЕНИЕ МАССИВА Y'
stop
endif
c-----------------------------------------------------------------
c проверка границы трубы с жидкостью
c-----------------------------------------------------------------
if ((IIY.ge.KR1).and.(IIY.le.KR2).and.(IIX.ge.KR3)) then
q(IIX,IIY)=200.
c-----------------------------------------------------------------
c проверка симметрии
c-----------------------------------------------------------------
elseif (((IIY.lt.KR1).or.(IIY.gt.KR2)).and.(IIX.eq.M)) then
q(IIX,IIY)=6
NN=NN+1
y(1,NN)=IIX
y(2,NN)=IIY
y(3,NN)=2*HX
y(4,NN)=KY
y(5,NN)=0
y(6,NN)=KY
c-----------------------------------------------------------------
c составление уравнений во внутренних точках фигуры
c-----------------------------------------------------------------
else
q(IIX,IIY)=5
NN=NN+1
y(1,NN)=IIX
y(2,NN)=IIY
y(3,NN)=HX
y(4,NN)=KY
y(5,NN)=HX
y(6,NN)=KY
endif
enddo
enddo
c-----------------------------------------------------------------
c присваивание начальных значений на границе фигуры
c------------------------------------------------------------------
KY=0
KR1=P/2+1
do IIY=1,P
if (IIY.le.KR1)then
q(1,IIY)=0
else
q(1,IIY)=500*KY-100
endif
KY=KY+K
enddo
do IIX=1,M
q(IIX,1)=0
q(IIX,P)=100
enddo
N=NN
end
c-----------------------------------------------------------------
c ПОДПРОГРАММА РЕШЕНИЯ СИСТЕМЫ УРАВНЕНИЙ МЕТОДОМ ЗЕЙДЕЛЯ
c
с integer N-входное количество уравнений
c real y(6,N)-входной массив уравнений,содержащий следующие поля:
c y(1,N)-номер точки по оси X
c y(2,N)-номер точки по оси Y
c y(3,N)-коэфициен уравнения для Q(y(1,N)-1,y(2,N))
c y(3,N)=h^2/(2*(h^2+k^2))
c y(4,N)-коэфициен уравнения для Q(y(1,N),y(2,N)-1)
c y(4,N)=k^2/(2*(h^2+k^2))
c y(5,N)-коэфициен уравнения для Q(y(1,N)+1,y(2,N))
c y(5,N)=h^2/(2*(h^2+k^2))
c y(6,N)-коэфициен уравнения для Q(y(1,N),y(2,N)+1)
c y(6,N)=k^2/(2*(h^2+k^2))
c integer M-число узлов по оси X
c integer P-число узлов по оси Y
c real Q(M,P)-входной массив начальных значений Y
c real Q(M,P)-выходной массив вычисленых значений Y
c real E-погрешность вычислений
c------------------------------------------------------------------
subroutine zeidel(N,y,M,P,q,E)
integer N,M,P,I,S
real y(6,N),q(M,P),E,EI,NEXTQ
c------------------------------------------------------------------
c вычисление коэфициента сходимости процесса
c mj=y(5,1)+y(6,1)
c и выражения
c km=mj/(1-mj)
C НО Т.К. MJ=0.5 ТО KM=1 И СЛЕДОВАТЕЛЬНО ЕГО