Температурный расчет с помощью вычислений информационной математики

То есть для
каждой строки
матрицы коэффициентов
системы
вычисляется
сумма модулей
коэффициентов,
лежащих правее
главной диагонали
:

и сумма модулей
коэффициентов,
лежащих левее
главной диагонали:

Для каждой
i-й
строки
(i =1,2,...,n
) вычисляется
отноше­ние

и в качестве
берется
максимальное
из этих отношений.
Чем меньше
окажется
,
тем большей
будет скорость
сходимости.

Для процесса
Зейделя справедлива
следующая
оценка погрешнос­ти
К-го приближения:

(i,j=1,2,...,n)

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

Если задаться
абсолютной
погрешностью
и
потребовать
выполнения
условия

(j=1,2,...,n)

то выполнится
и условие

(i=1,2,3,...,n),

то есть заданная
степень точности
на К-й итерации
будет достигнута.
На практике
это означает,
что после каждой
итерации необходимо
вы­делить тот
корень, изменение
которого по
сравнению с
предыдущим
значением
оказалось
наибольшим
по модулю. Модуль
приращения
этого корня
необходимо
умножить на
и сравнить
результат с
выбран­ной
абсолютной
погрешностью.

Достаточные
условия сходимости
процесса Зейделя

Если модули
коэффициентов
системы
удовлетворяют
хотя бы одному
из условий

(i,j=1,2,3,...,n)

то процесс
Зейделя для
соответствующей
приведенной
системы сходит­ся
к её единственному
решению при
любом выборе
начального
вектсра y(ср)
Такие системы
называют системами
с диагональным
преоблада­нием.

Метод Зейдедя
имеет свойство,
позволяющее
обеспечить
сходимость
процесса для
любых систем
уравнений с
неособенной
матрицей
коэфициентов.

Если обе
части систем
с неособенной
матрицей коэфициентов
А=[aij]
умножить слева
на транспонировнную
матриц A*[aij]
, то будет
получена новая,
равносильная
исходной система,
которая называется
нормальной.
Процесс Зейделя
для приведенной
системы, полученной
из нормальной,
всегда сходится
независимо
от выбора нача
льного приближения.

Блок
схема.

Подпрограмма
метода Зейделя.

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 И
СЛЕДОВАТЕЛЬНО
ЕГО МОЖНО ОПУСТИТЬ

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

c------------------------------------------------------------------

c
ПОДПРОГРАММА
АЛФАВИТНО-ЦИФРОВОГО,МОЗАИЧНОГО

c
ВЫВОДА
РЕЗУЛЬТАТА

c
integer M-число узлов
по оси X

c
integer P-число узлов
по оси Y

c
real Q(M,P)-входной
массив значений
Y

c

c

c------------------------------------------------------------------

subroutine outdata(M,P,q)

character
a(11)/'.','+','*','','','-','-','-','','-','-'/

integer M,P,I,J

real
q(M,P)

do
J=P,1,-1

print '(400A2)',(a(int(q(I,J)/21)+1),I=1,M),

+
(a(int(q(I,J)/21)+1),I=M-1,1,-1)

enddo

do
I=1,10

print *,'''',a(I),'''','---> от
',20*(I-1),', до ',

+
20*I,'(включительно)'

enddo

end

c------------------------------------------------------------------

c
ПОДПРОГРАММА
ВЫЧИСЛЕНИЯ
ОШИБКИ

c
real q-массив значений
Y с шагом =2*h

c
real qq-массив значений
Y с шагом =h

c
real E-значение
погрешности

c

c------------------------------------------------------------------

subroutine mistake(M,P,q,qq,E)

integer M,P,iq,jq,iqq,jqq

real
qq(M,P),q(int(M/2)+1,int(P/2)+1),max,E,other

max=0

iq=0

do
iqq=1,P,2

iq=iq+1

jq=0

do
jqq=1,M,2

jq=jq+1

other=abs(q(jq,iq)-qq(jqq,iqq))

if
(other.gt.max)max=other

enddo

enddo

print *,M,' ',P,' ',max/3

if
(max/3.lt.E) then

call outdata(M,P,qq)

Stop

endif

end

c------------------------------------------------------------------

c
ОСНОВНАЯ
ПРОГРАММА

c

c

c------------------------------------------------------------------

integer N/90000/,M,P,flag/0/

real
y(6,90000),q(300,300),H/.05/,K/.05/,E/.5/,qq(300,300)

real
EZ/.01/

c
print *,'Введите шаг
вдоль оси X '

c
read (*,*)H

c
print *,'Введите шаг
вдоль оси Y '

c
read (*,*)K

c
print *,'Введите
точность вычислений
'

c
read (*,*)E

M=.2/H+1

P=.4/K+1

call
mkr(H,K,N,y,M,P,q)

call
zeidel(N,y,M,P,q,EZ)

111
H=H/2

K=K/2

M=.2/H+1

P=.4/K+1

N=90000

if
(flag.eq.0)then

flag=1

call mkr(H,K,N,y,M,P,qq)

call zeidel(N,y,M,P,qq,EZ)

call mistake(M,P,q,qq,E)

else

flag=0

call mkr(H,K,N,y,M,P,q)

call zeidel(N,y,M,P,q,EZ)

call mistake(M,P,qq,q,E)

endif

goto
111

end

Литература.

1.
И.С.Березин,Н.П.Жидков
’Методы
вычислений’,том
1,М.,1966,632 стр.

2.’
Численные
методы решения
задач на ЭВМ
’ , Учебное
пособие , Г.Н.Рубальченко
, К. ,
1989
, 148 стр.

3.’Справочник
языка ФОРТРАН’
, М.,1996 ,106 стр.

Температурный
расчет с помощью
вычислений
информационной
математики.