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

Постановказадачи.

Подлинной квадратногосечения трубетечет горячаяжидкость. Трубанаполовинупогружена вледяную ванну,так, что температуранижней половиныповерхноститрубы равна 0⁰С. Верхняя плоскостьтрубы имеетпостояннуютемпературу100 ⁰С. На участкемежду ледянойванной и верхнейплоскостьютемпературанаружной поверхноститрубы изменяетсялинейно повысоте от 0⁰ С до100 ⁰С. Жидкостьвнутри трубыимеет температуру200 ⁰С.

Рис. 3.

Распределениетемпературывтеле трубыудовлетворяетуравнению

Спогрешностьюне более 0,5⁰ С вычислитьраспределениетемпературыв теле трубы.

Математическаяформулировказадачи.

Решитьдиф.уравнениев частныхпроизводных:

сзадаными началинымиусловиями награницах областидифференцирования.

Прирешении уравненияприблизительнозаменю производныевторого порядкаконечно-разностнымиотношениями:

врезультатечего диф.уравнениепреобразуетсяв 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 стр.

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