Постановказадачи.
Подлинной квадратногосечения трубетечет горячаяжидкость. Трубанаполовинупогружена вледяную ванну,так, что температуранижней половиныповерхноститрубы равна 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 и ограничившисьчетырьмя членами:
где
Аналогичнонайдем значениеф-ции и в точке
где
Подставляяполучим:
Такимобразом,производная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и т.д.
Тогда получимэквивалентнуюсистему
где
и
Такую системубудем вдальнейшемназывать приведенной.
Метод Зейделязаключаетсяв следующем.Выбрав векторначальногоприближения
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 -го корнясистемы в К-мприближенииот точногозначения тогоже корня
Если задатьсяабсолютнойпогрешностью
то выполнитсяи условие
то есть заданнаястепень точностина К-й итерациибудет достигнута.На практикеэто означает,что после каждойитерации необходимовыделить тоткорень, изменениекоторого посравнению спредыдущимзначениемоказалосьнаибольшимпо модулю. Модульприращенияэтого корнянеобходимоумножить на
Достаточныеусловия сходимостипроцесса Зейделя
Если модуликоэффициентовсистемы удовлетворяютхотя бы одномуиз условий
то процессЗейделя длясоответствующейприведеннойсистемы сходитсяк её единственномурешению прилюбом выбореначальноговектсра 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 стр.
Температурныйрасчет с помощьювычисленийинформационнойматематики.