Температурный расчет с помощью вычислений информационной математики
Постановказадачи.
Подлинной квадратногосечения трубетечет горячаяжидкость. Трубанаполовинупогружена вледяную ванну,так, что температуранижней половиныповерхноститрубы равна 00С. Верхняя плоскостьтрубы имеетпостояннуютемпературу100 0С. На участкемежду ледянойванной и верхнейплоскостьютемпературанаружной поверхноститрубы изменяетсялинейно повысоте от 00 С до100 0С. Жидкостьвнутри трубыимеет температуру200 0С.
Рис. 3.
Распределениетемпературы
втеле трубыудовлетворяетуравнению
Спогрешностьюне более 0,50 С вычислитьраспределениетемпературыв теле трубы.
Дискретизация | Методконечныхразностей | + |
задачи | Методконечныхэлементов | |
Решение | МетодГаусса | |
системы | МетодЗейделя | + |
линейных | Методпоследовательнойверхней релаксации | |
уравнений | Методрелаксацияпо строкам | |
Вывод | Библиотечнаяграфическаяподпрограмма | |
результатов | Алфавитно-цифровой,мозаичный | + |
Математическаяформулировказадачи.
Решитьдиф.уравнениев частныхпроизводных:
сзадаными началинымиусловиями награницах областидифференцирования.
Прирешении уравненияприблизительнозаменю производныевторого порядкаконечно-разностнымиотношениями:
врезультатечего диф.уравнениепреобразуетсяв 5-ти диаганальнуюсистему алгеброическихуравнений n-гопорядка.
Системуалгеброическихуравнений будурешать методомЗейделя.
Погрешностьрешения задачинайду по формуле:
где,
и 
-решения,полученныедля одной и тойже точки с разнымишагами.Функциональнаясхема.
Метод конечныхразностей.
Описаниеметода.
Так названметод решениякраевых задач,основанныйна приближеннойзамене производных,входящих вдифференциальныеуравненияи краевые условия,нонечно-разностнымиотношениями.Эта заменапозволяетсвести краевуюзадачу к задачерешения системыалгебраическихуравнений.
Конечныеразности ипроизводные.Пустьнекотораяфункция y(x)задана на отрезке[a,b].Будем считать,что она непрерывнаи многократнодифференцируемана этом отрезке.Разделим отрезокна равные частидлиною hи обозначимточки деленияx0,x1,...,xi,...,xn.Значенияфункции в этихточках обозначимсоответственноy0,y1,...,yi,...,yn.Первойцентральнойразностью вi-йточке(i=1,2,...,n-1) называютразность:
С помощьюэтой разностиможно приближенновычислитьзначение первойпроизводнойу`в i-йточке.
Разложимфункцию y(x)в степеннойряд. приняв зацентр разложенияточкуxi и ограничившисьчетырьмя членами:
где
Аналогичнонайдем значениеф-ции и в точке
,отстоящейот центра разложенияна шаг (-h):
где
.Подставляяполучим:
Такимобразом,производнаяy`приближоннозаменяетсяконечно-разностнымотношениемс ошибкой порядкаh*h:
Второйцентральнойразностью ф-цииy(x)в i-йточке называютвеличину:
Спомощью этойразности можноприближонновычислитьзначение второйпроизводнойy``в i-йточке.Используемтеперь 5 членовразложенияв ряд Тейлора:
Такимобразом,втораяпроизводнаяy``с ошибкой порядкаh*hможет бытьприближоннозамененаконечно-разностнымотношением:
При определенииразностей вi-и точке использовалисьзначения функциив точках, расположенныхсимметричноотносительноxi. Поэтомуэти разностиназываютсяцентральными.
Существуюттакже левыеи правые разности,использующиеточки, расположенныесоответственнолевее и правееточки xi.С помощью этихразностей можнотакже приближенновычислятьзначения производных,но погрешностьпри этом будетбольше-порядкаh.
Разностныесистемы уравненийсоставляютсяв следующемпорядке.
1. Исходноедифференциальноеуравнениепреобразуютк такой форме,чтобы затемполучить изнего наиболеепростую разностнуюсистему уравнений.При этом учитывают,что коэффициентыпри производныхвойдут в разностнуюсхему одновременнов несколькоее членов изатем будутраспространенына всю системууравнений.Поэтому желательноиметь единичныекоэффициентыпри производныхв исходномуравнении.
2. Наинтервалеинтегрированияисходногоуравненияустанавливаютравномернуюсетку с шагомh и записываютразностнуюсхему, приближеннозаменяя производныесоответствующимицентральнымиконечно-разностнымиотношениями.
3.Применяяразностнуюсхему для узловсетки записываютразностныеуравнения. Приэтом можно получить уравнениясодержащиетак называемыевнеконтурныенеизвестные,то есть неизвестныев точках, лежащихза пределамиустановленнойсетки.
4.В разностнойформе записываюткраевые условияи составляютполную системуразностныхуравнений.
Оценкапогрешностирешения краевойзадачи
Решениеразностнойсистемы уравненийдает приближенноерешениекраевой задачи.Поэтому возникаетвопрос о точностиэтого приближенногорешения.
Для линейныхкраевых задачдоказана теоремао том, что порядокточности решениякраевой задачине ниже порядкаточностиаппроксимациипроизводныхконечно-разностнымиотношениями.Оценку погрешностипроизводятприемом Рунге.Краевую задачурешают дважды:с шагом сеткиh ис шагом сеткиH=kh, погрешностьрешения с малымшагом hоцениваютпо формуле:
где y(h)и y(H)- решения,полученныедля одной и тойже точки-xi отрезкаинтегрированияс разными шагами.ОтносительнуюпогрешностьEоценивают поформуле:
Если присоставленииразностнойсистемы уравненийиспользуютсялевые или правыеразности, топогрешностьрешения будетвыше, порядка0(h),и для ее оценкив формулахследует заменитьk*k на k.
Применениеметода конечныхразностей длярешения уравненийв частных проиэводных
Дляпримененияразностногометода в областиизменениянезависимыхпеременныхвводят некоторуюсетку. Всепроизводные,входящие вуравнение икраевые условия,заменяют разностямизначений функциив узлах сеткии получаюттаким образомалгебраическугосистему уравнений.Решая эту систему,находят приближенноерешение задачив узлах сетки.
Блоксхема.
ПодпрограммаМКР.
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
realy(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
doIIY=2,P-1
doIIX=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
doIIY=1,P
if(IIY.le.KR1)then
q(1,IIY)=0
else
q(1,IIY)=500*KY-100
endif
KY=KY+K
enddo
doIIX=1,M
q(IIX,1)=0
q(IIX,P)=100
enddo
N=NN
end
ТЕСТ
Длятестированиясоставлю разностнуюсистему с шагомвдоли оси Xи Y=0.05
Неизвестныезначения вузлах матрицынаходящихсявнутри фигурывысчитываютсяпо формуле:
Неизвестныезначения вузлах матрицынаходящихсяна оси симметриивысчитываютсяпо формуле:
МЕТОД ЗЕЙДЕЛЯ.
Метод Зейделяотносится кчислу итерационныхметодов, в которыхпринципиальноотсутствуетфактор накопленияпогрешностей.Поэтому оншироко применяетсядля решениябольших системуравнений.Будем рассматриватькорни решаемойсистемы каккомпонентынекотороговектора у. Основнаяидея всехитерационныхметодов заключаетсяв том, что беретсяприближенноезначение векторау и по формулам,составленнымна основаниирешаемых уравнений,вычисляетсяновое приближенноезначение векторау .Назовем этиприближенныезначенияy(k) иy(k+1) соответственно.Посколькуисходноеприближениевыбиралосьпроизвольно,то у(k+1)всвою очередьможет послужитьисходным дляполучения потем же формуламнового приближенияy(k+2). Очевидно,этот процессможно продолжатьсколь угоднодолго. Говорят,что процесситераций сходится,если получаемаяпри этом последовательностьвекторов у(k) (к=0,1,2,...}имеет своимпределом векторy,являющийсяточным решениемсистемы:
На практикеневозможнодостигнутьэтого предела,но можно приблизитьсян нему с любойнаперед заданнойточностью. Таки поступаютзадаются некоторойпогрешностью,вектором начальногоприближенияи получаютпоследовательныеприближениядо тех пор, покадействительнаяпогрешностькорней не станетменьше заданной.
Различныеметоды отличаютсядруг от другаспособом вычисленияочередногоприближения,но во всех методахсуществуютдве главныепроблемы:
обеспечениесходимостипроцесса итераций;
оценкадостигнутойпогрешности.
Пусть даналинейная система
Предполагая,что диагональныекоэффициенты
aii0(i=1,2,..,n)
разрешимпервое уравнениеотносительноy1 ,второе- относительноy2и т.д.
Тогда получимэквивалентнуюсистему
где
при ijи
при i=j (i,j=1,2,...,n)Такую системубудем вдальнейшемназывать приведенной.
Метод Зейделязаключаетсяв следующем.Выбрав векторначальногоприближения
y(ср)=(y1,y2,...,yn)
подставимего компонентыв правую частьпервого уравнениясистемыи вычислимпервую компоненту y`1нового вектораy`(ср) . В правуючасть второгоуравненияподставимкомпоненты(y`1,y2,y3,...,yn)и вычислимвторую компонентуy`2'новоговектора. В третьеуравнениеподставим(y`1,y`2,y3,...,yn) и т.д.Очевидно,подстановкойв каждоеуравнение мы,дойдя до последнегоуравнения,обновим всекомпонентыисходноговектора и получимпервое приближениек решению
y`(ср)=(y`1,y`2,y`3,...,y`n)
Далее, взяв вкачестве исходноговектор у`(ср) , выполнимвторуюитерациюи получим y``(ср).Этот процессбудем продолжатьдо достижениязаданной степениточности.
Оценкапогрешностиприближенийпроцесса Зейделя
Для оценкипогрешностипрежде всеговычисляютпоказательскорости сходимости
То есть длякаждой строкиматрицы коэффициентовсистемы
вычисляетсясумма модулейкоэффициентов,лежащих правееглавной диагонали 
:
и сумма модулейкоэффициентов,лежащих левееглавной диагонали:
Для каждой 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
realy(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
doI=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
realy(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
doIIY=2,P-1
doIIX=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
doIIY=1,P
if(IIY.le.KR1)then
q(1,IIY)=0
else
q(1,IIY)=500*KY-100
endif
KY=KY+K
enddo
doIIX=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
realy(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
doI=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)
charactera(11)/'.','+','*','','','-','-','-',''
,'-','-'/integer M,P,I,J
realq(M,P)
doJ=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
doI=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
realqq(M,P),q(int(M/2)+1,int(P/2)+1),max,E,other
max=0
iq=0
doiqq=1,P,2
iq=iq+1
jq=0
dojqq=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/
realy(6,90000),q(300,300),H/.05/,K/.05/,E/.5/,qq(300,300)
realEZ/.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
callmkr(H,K,N,y,M,P,q)
callzeidel(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
goto111
end
Литература.
1.И.С.Березин,Н.П.Жидков’Методывычислений’,том1,М.,1966,632 стр.
2.’Численныеметоды решениязадач на ЭВМ’ , Учебноепособие , Г.Н.Рубальченко, К. ,1989, 148 стр.
3.’Справочникязыка ФОРТРАН’, М.,1996 ,106 стр.
Температурныйрасчет с помощьювычисленийинформационнойматематики.