Многомерные нестационарные задачи теплопроводности. Примеры реализации разностного решения в среде Х

Государственное образовательное учреждение высшего профессионального образования «Московский государственный технический университет имени Н.Э. Баумана»

Калужский филиал

Факультет «Фундаментальных Наук»

Кафедра «Программного Обеспечения ЭВМ, Информационных Технологий и Прикладной Математики»

РАСЧЕТНО-ПОЯСНИТЕЛЬНАЯ ЗАПИСКА К КУРСОВОЙ РАБОТЕ

Тема: «Многомерные нестационарные задачи теплопроводности. Примеры реализации разностного решения в среде Х»

Студент: Свободин А.М. (_______)

Группа ИТД — 82

Преподаватель: Гинзгеймер С.А. (_______)

2010 г

ОГЛАВЛЕНИЕ

• 1. Техническое задание
• 1.1 Введение
• 1.2 Основания для разработки
• 1.3 Назначение разработки
• 1.4 Требования к программному проекту
• 1.4.1 Требования к функциональным характеристикам
• 1.4.2 Требования к надёжности
• 1.4.3 Условия эксплуатации
• 1.4.4 Требования к составу и параметрам технических средств
• 1.4.5 Требования к информационной и программной совместимости
• 1.5 Требования к программной документации
• 1.6 Стадии и этапы разработки
• 2. Исследовательская часть
• 2.1 Постановка задачи
• 2.2 Выбор средства разработки
• 2.3 Общие сведения
• 2.3.1 Температурное поле
• 2.3.2 Разностные схемы
• 2.2.3.2 Аппроксимация и порядок аппроксимации
• 2.3.2.3 Устойчивость
• 2.3.2.4 Явные схемы
• 2.3.2.5 Неявные схемы
• 3. Конструкторская часть
• 3.1 Общие сведения
• 3.2 Составление численного решения
• 3.2.1 Метод переменных направлений
• 3.2.2 Метод дробных шагов
• 3.4 Назначение программы
• 3.5 Условия выполнения программы
• 3.6 Тестирование программы
• 3.7 Обращение к программе
• 3.8 Структура программы
• 4. Технологическая часть
• 4.1 Входные данные
• 4.2 Выходные данные
• 4.3 Основные используемые функции
• 4.4 Результат работы скрипта
• 4.5 Анализ
• Литература
• Приложение 1

1. ТЕХНИЧЕСКОЕ ЗАДАНИЕ

1.1 Введение

программа аппроксимация температурный дробный

Данная курсовая работа должна представлять собой набор скриптов пакета Maple, предназначенный для расчета температурного поля плоской прямоугольной пластины.

1.2 Основания для разработки

Данное приложение должно быть разработано на основании курсовой работы по курсу «Математическое моделирование» на тему «Многомерные нестационарные задачи теплопроводности».

1.3 Назначение разработки

Приложение должно предназначаться для отображения температурного поля плоской прямоугольной пластины с возможностями изменения характеристики пластины, начальных условий и проведения анализа процесса.

1.4 Требования к программному проекту

1.4.1 Требования к функциональным характеристикам

Приложение должно содержать и выполнять следующие действия:

· возможность ввода и изменения характеристик пластины и начальных условий процесса;

· расчет температурного поля плоской прямоугольной пластины;

· вывод графиков температурного поля плоской пластины.

1.4.2 Требования к надёжности

Требования к надёжности обеспечиваются использованием лицензионного стандартного программного обеспечения.

1.4.3 Условия эксплуатации

Использование скрипт-программы может осуществлять пользователь, владеющий основными навыками работы с персональным компьютером.

1.4.4 Требования к составу и параметрам технических средств

Для нормального функционирования приложения необходимы:

· ОС Windows 9x/ME/2000/XP

· Процессор i80486

· 64Мб RAM

· Не менее 3 Мб свободного места на жестком диске

· Пакет Maple 13

1.4.5 Требования к информационной и программной совместимости

Среда проектирования, используемая при разработке данного проекта — Maple 13. Для работы с приложением требуется операционная система Windows 9x/ME/2000/XP.

1.5 Требования к программной документации

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

1.6 Стадии и этапы разработки

В стадию разработки технического задания включаются следующие этапы:

№ этапа	Наименование этапа	Сроки проведения	Содержание работ
1.	Обоснование необходимости разработки программы	08.02.2010 -16.02.2010	Проведение постановки задачи.
2.	Научно-исследовательские работы	17.02.2010 — 13.03.2010	Определение структуры входных и выходных данных, выбор методов решения задач, определение требований к техническим средствам.
3.	Разработка и утверждение технического задания	14.03.2010 — 28.03.2010	Определение требований к программе, установка сроков её разработки.

В стадию разработки технического проекта входят этапы:

№ этапа	Наименование этапа	Сроки проведения	Содержание работ
1.	Разработка технического проекта	29.03.2010 — 20.04.2010	Разработка алгоритма решения задачи, уточнение структур входных и выходных данных.
3.	Утверждение технического проекта	20.04.2010 — 10.05.2010	Разработка пояснительной записки, согласование и утверждения технического проекта.

Исполнителем программного проекта является Свободин А.М.

Данный проект должен быть утверждён до 15.05.2010 г.

2. Исследовательская часть

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

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

Нестационарное уравнение теплопроводности плоской прямоугольной пластины:

(1)

— функция, возвращающая температуру определенной точки пластины в заданный момент времени.

— коэффициент теплопроводности вдоль оси х.

— коэффициент теплопроводности вдоль оси у.

— функция источника.

Программа должна позволять изменять геометрический размер пластины, коэффициенты теплопроводности, а также начальные условия.

Среда проектирования — Maple 13.

2.2 Выбор средства разработки

Для реализации курсовой работы была использована система Maple 13.

Maple — cсистема компьютерной алгебры, ориентированная на сложные математические вычисления, визуализацию данных и моделирование. Система Maple предназначена для символьных вычислений, хотя имеет ряд средств и для численного решения дифференциальных уравнений и нахождения интегралов. Обладает развитыми графическими средствами. Имеет собственный язык программирования, частично подобный Pascal.

Начинающим пользователям система Maple 13 предоставляет средства обучения и инструменты для пошагового решения наиболее сложных задач. Специалистам предоставляется удобный доступ к дополнительным ресурсам для быстрого получения ответов на различные вопросы. Продукт Maple 13 интегрирован с ведущими САПР-системами, включая популярную платформу NX. Maple 13 предлагает усовершенствованные и высокопроизводительные средства численных и символьных вычислений с возможностью представления математических результатов в виде подробной технической документации. Система обеспечивает высокое качество графики и анимации и предоставляет комплекс инструментов для редактирования документов и средства управления визуализацией результатов

2.3 Общие сведения

2.3.1 Температурное поле

Рассмотрим однородные тела, т.е. такие тела, которые обладают одинаковыми физическими свойствами по всем направлениям. При передачи теплоты в твердом теле, температура тела будет изменяться по всему объему тела и во времени. Совокупность значений температуры в данный момент времени для всех точек изучаемого пространства называется температурным полем:

,

где:

T -температура тела;

x, y, z — координаты точки;

— время.

Такое температурное поле называется нестационарным.

Если температура тела функция только координат и не изменяется с течением времени, то температурное поле называется стационарным:

, .

2.3.2 Разностные схемы

2.3.2.1 Общие сведения

Разностная схема — разностный метод приближенного решения какого-либо дифференциального уравнения в частных производных или применения дифференциального оператора. Разностные схемы применяются к функциям, заданным на какой-либо сетке.

Разностная схема, как правило, использует уравнения, связывающие несколько соседних точек результата и исходных данных (результата на предыдущих шагах в случае дифференциального уравнения). Решение этих уравнений позволяет найти приближенное решение.

2.2.3.2 Аппроксимация и порядок аппроксимации

Запишем дифференциальную задачу в операторной форме

LU = f, где L — один из дифференциальных операторов

U(x,y) — искомая функция удовлетворяющая дифференциальной задаче; f — входные данные (т.е. начальные и краевые условия, правые части и т.п.). Операторная форма описывает дифференциальную задачу в узлах сетки, а операторная форма описывает конечно-разностную схему на точном решении U(x,t),т.е. в конечно-разностной схеме вместо сеточных значений сеточной функции подставлены точные (неизвестные) значения искомой функции. Операторная форма конечно-разностной схемы имеет вид .

Введём норму сеточной функции, например, с помощью выражения

Определение. Конечно-разностная схема аппроксимирует дифференциальную задачу на точном решении, если какая-либо норма разности ? стремится к нулю при

Определение. Конечно-разностная схема аппроксимирует дифференциальную задачу на точном решении с порядком p по времени и порядком q по пространственной переменной, если какая-либо норма разности удовлетворяет равенству

Таким образом, если конечно-разностная схема аппроксимирует дифференциальную задачу, то речь идёт о близости дифференциального и конечно-разностного операторов в узлах сетки.

2.3.2.3 Устойчивость

Пусть в конечно-разностной схеме входные данные получили возмущения и стали . Тогда сеточная функция также получит возмущение и станет .

Определение. Конечно-разностная схема устойчива по входным данным, если найдется такая ограниченная константа K, не зависящая от сеточных характеристик и входных данных, что выполняется неравенство (4)

Таким образом, понятие устойчивости интерпретируется следующим образом: конечно-разностная схема устойчива, если для малых возмущений входных данных (начально-краевых условий и правых частей) конечно-разносная схема обеспечивает малые возмущения сеточной функции , т.е. решение с помощью конечно-разностной схемы находится под контролем входных данных.

Если во входные данные входят только начальные условия или только краевые условия, или только правые части, то говорят об устойчивости соответственно по начальным условиям, по краевым условиям или по правым частям

Определение. Конечно-разностная схема абсолютно (безусловно) устойчива, если неравенство (4) выполняется при любых значениях сеточных характеристик ф и h, т.е. на шаги сетки не накладывается никаких ограничений.

Определение. Конечно-разностная схема условно устойчива, если неравенство (4) выполняется для сеточных характеристик ф и h, на которые накладываются определённые ограничения.

2.3.2.4 Явные схемы

Рассмотрим первую начально-краевую задачу для волнового уравнения. На пространственно-временной сетке будем аппроксимировать дифференциальное уравнение одной из следующих конечно-разностных схем:

с шаблоном на рисунке 1 и

с шаблоном на рисунке 2

рис 1.

рис 2.

При этом схема (1) является явной. С ее помощью решение , определяется сразу, поскольку значения сеточных функции на нижних временных слоях должны быть известны. В соответствии с шаблоном для этой схемы порядок аппроксимации равен двум, как по пространственной, так и по временной переменной. При этом явная конечно-разностная схема (1) для волнового уравнения условно устойчива с условием , накладываемым на сеточные характеристики ф и h.

2.3.2.5 Неявные схемы

Схема (2) является неявной схемой и обладает абсолютной устойчивостью. Ее можно свести к СЛАУ с трехдиагональной матрицей, решаемой методом прогонки.

Для определенияможно воспользоваться простейшей аппроксимацией второго начального условия.

Откуда для искомых значений получаем следующее выражение:

Неявные схемы обычно являются устойчивыми.

3. КОНСТРУКТОРСКАЯ ЧАСТЬ

3.1 Общие сведения

Данная курсовая работа представляет собой набор скриптов, созданный в среде Maple 13, с полным её описанием.

3.2 Составление численного решения

Уравнение теплопроводности пластины будет иметь вид:

с начальными и граничными условиями:

При численном решении многомерных задач математической физики исключительно важным является вопрос об экономичности используемых методов.

Конечно — разностную схему будем называть экономичной, если число выполняемых операций (операций типа умножения) пропорционально числу узлов сетки.

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

Из экономичных конечно-разностных схем, получивших наибольшее распространение, является схема метода переменных направлений и схема метода дробных шагов. Все эти методы обычно называются общим термином — методы расщепления.

3.2.1 Метод переменных направлений

В схеме метода переменных направлений (МПН), как и во всех методах расщепления, шаг по времени ф разбивается на число независимых пространственных переменных (в двумерном случае — на два ). На каждом дробном временном слое один из пространственных дифференциальных операторов аппроксимируется неявно (по соответствующему координатному направлению осуществляются скалярные прогонки ), а остальные явно. На следующем дробном шаге следующий по порядку дифференциальный оператор аппроксимируется неявно, а остальные — явно и т.д. В двумерном случае схема метода переменных направлений имеет вид

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

3.2.2 Метод дробных шагов

В отличие от МПН метод дробных шагов (МДШ) использует только неявные конечно-разностные операторы, что делает его абсолютно устойчивым в задачах, не содержащих смешанные производные. Он обладает довольно значительным запасом устойчивости и в задачах со смешанными производными.

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

К недостаткам МДШ относятся следующие: на каждом дробном шаге достигается частичная аппроксимация, полная аппроксимация достигается на последнем дробном шаге, т.е. имеет место суммарная аппроксимация; схема имеет первый порядок точности по времени.

3.3 Назначение программы

Данный набор скриптов создан в среде Maple 13. Он позволяет пользователю получать графическое изображение температурного поля прямоугольной пластины, рассчитанного разными методами и с помощью возможности изменять характеристики пластины и начальных условий проводить анализ.

3.4 Условия выполнения программы

Программа может быть запущена на любом компьютере, на котором установлена одна из операционных систем Windows 9x/ME/NT/2000/XP. Набор скриптов занимает 90 КБ дискового пространства. Для их выполнения достаточно 64 Мб оперативной памяти.

3.5 Тестирование программы

Для выполнения тестирования решения следует произвести следующие действия:

1. Выбрав требуемый метод решения запустить файл «Дробных шагов.mws» либо «Переменных направлений.mws».

2. Задать параметры пластины и начальные условия процесса.

3. Просмотреть графики температурного поля пластины.

4. Выполнить выход из Мaple.

Тестирование скрипта окончено.

3.6 Обращение к программе

Для обращения к скрипту необходимо запустить пакет Maple 13 открыть файл «Дробных шагов.mws» либо «Переменных направлений.mws».

3.7 Структура программы

Приложение состоит из трех файла:

1. Дробных шагов.mws,

2. Переменных направлений.mws,

3. diff.mws;

соответственно выполняющих расчеты с применением метода Дробных шагов и Переменных направлений и производящих сравнение результатов расчета этих двух методов.

4. ТЕХНОЛОГИЧЕСКАЯ ЧАСТЬ

4.1 Входные данные

Входными данными являются следующие характеристики пластины:

lx — длина пластины;

ly — ширина пластины;

tim — время в течении которого наблюдается процесс теплопроводности;

Nx, Ny, Nt — размерность сетки разностного решения;

f — функция задающая температуру источника тепла;

ksi — функция задающая начальную температуру;

fi1, fi2, fi3, fi4 — функции, задающие температуру на границе пластины;

a — квадрат скорости распространения тепла вдоль оси х;

b — квадрат скорости распространения тепла вдоль оси у;

4.2 Выходные данные

Выходными данными являются графики температурного поля пластины.

4.3 Основные используемые функции

LinearSolve(M) — функция, решающая алгебраическую систему уравнений (матрица М) surfdata(g, r1, r2, options) — создает 3D поверхность используя входные данные g- матрица значений, r1,r2- диапазон выбираемых значений.

4.4 Результат работы приложения

lx=2*Pi;

ly=2*Pi;

a=1;

b=1;

f:=(x,y,t)->-x*y*sin(mu*t):

ksi:=(x,y)->x*y:

fi1:=(y,t)->0.0001:

fi2:=(y,t)->y*cos(mu*t):

fi3:=(x,t)->0.0001:

fi4:=(x,t)->x*cos(mu*t):

tim=Pi/60;

Рис 3. Аналитическое решение

Рис 4. Численное решение (схема переменных направлений)

Рис 5. Разность между аналитическим решением и численным (переменных направлений)

Рис 6. Численное решение (схема дробных шагов)

Рис 7. Разность между аналитическим решением и численным (дробных шагов)

Рис 8. Разность между численными решениями

tim=Pi/4;

Рис 9. Аналитическое решение

Рис 10. Численное решение (схема переменных направлений)

Рис 11. Разность между аналитическим решением и численным (переменных направлений)

Рис 12. Численное решение (схема дробных шагов)

Рис 13. Разность между аналитическим решением и численным (дробных шагов)

Рис 14. Разность между численными решениями

4.5 Анализ

В ходе анализа получившихся решений (рис 3 — рис 14) было обнаружено что численное и аналитические решения совпадают при малых t, с ростом же t увеличивается ошибка на границе пластины. Аналогично дело обстоит с разницей между численными решениями. Отсюда можно сделать вывод об удовлетворительной устойчивости схем расщепления для первоначальной оценки процессов теплопередачи для коротких промежутков времени.

ЛИТЕРАТУРА

1) Лыков А.В. Теория теплопроводности

2) Полянин А.Д. Справочник по линейным уравнениям математической физики

3) Калиткин Н.Н. Численные методы

4) Самарский А.А. Введение в численные методы

5) www.wikipedia.org

6) www.intuit.ru

ПРИЛОЖЕНИЕ 1

Файл «Переменных направлений.mws»

#u_t = a*u_xx + b*u_yy + f(x,y,t),

#u(0,y,t) = fi_1(y,t)=0,

#u(lx,y,t) = fi_2(y,t)=y*cos(mu*t),

#u(x,0,t) = fi_3(x,t)=0,

#u(x,ly,t) = fi_4(x,t)=x*cos(mu*t),

#u(x,y,0) = ksi(x,y)=0,

#f(x,y,t) = x*y*sin(mu*t)).

restart:

#Коэффициенты теплопроводности#

a:=1:

b:=1:

mu:=1:

Nx:=20:

Ny:=20:

Nt:=30:

lx:=2*Pi:

ly:=2*Pi:

tim:=Pi/2:

fi1:=(y,t)-0.0001:

fi2:=(y,t)-y*cos(mu*t):

fi3:=(x,t)-0.0001:

fi4:=(x,t)-x*cos(mu*t):

ksi:=(x,y)-x*y:

f:=(x,y,t)—x*y*sin(mu*t):

hx:=lx/Nx:

hy:=ly/Ny:

tau:=tim/Nt:

with (LinearAlgebra):

#Задаем кубическую матрицу решения#

U:=vector(Nt+1):

for k from 1 to Nt+1 do

U[k]:=Matrix(Nx+1,Ny+1,0);

end do:

# Заполняем граничные условия

for k from 2 to Nt+1 do

for i from 1 to Nx+1 do

U[k][i, 1]:=fi3((i-1)*hx, (k-1)*tau);

U[k][i, Ny+1]:=fi4((i-1)*hx, (k-1)*tau);

end do:

for j from 1 to Ny+1 do

U[k][1, j]:=fi1((j-1)*hy, (k-1)*tau);

U[k][Nx+1, j]:=fi2((j-1)*hy, (k-1)*tau);

end do:

end do:

# Закончили заполнять граничные условия

# Заполняем начальные условия

for i from 1 to Nx+1 do

for j from 1 to Ny+1 do

U[1][i, j]:=ksi((i-1)*hx, (j-1)*hy);

end do:

end do:

# Закончили заполнять начальные условия

# Начинаем решать-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|

U12:=Matrix(Nx+1,Ny+1):

for k from 1 to Nt do

# Шаг 1 решаем ряды х, поднимаясь по у

====================================================

A_const:=a/(hx^2);

B_const:=-2*a/(hx^2)-2/tau;

C_const:=a/(hx^2);

for j from 2 to Ny do

SST:=Matrix(Nx+1, Nx+2,0);

for i from 2 to Nx do

SST[i,i-1]:=A_const;

SST[i,i]:=B_const;

SST[i,i+1]:=C_const;

SST[i, Nx+2]:=-b*(U[k][i, j+1]-2*U[k][i, j]+U[k][i, j-1])/(hy^2) — f((i-1)*hx,(j-1)*hy,((k-1)+1/2)*tau) -2*U[k][i, j]/tau;

end do:

SST[1,1]:=C_const;

SST[Nx+1, Nx+1]:=A_const;

SST[1, Nx+2]:=U[k][1,j];

SST[Nx+1, Nx+2]:=U[k][Nx+1,j];

xsol:=LinearSolve(SST):

for i from 2 to Nx do

U12[i,j]:=xsol[i]:

end do:

# Закончили решаеть вдоль у

# Заполняем новые граничные условия

U12[1, j]:=fi1((j-1)*hy,((k-1)+1/2)*tau);

U12[Nx+1, j]:=fi2((j-1)*hy, ((k-1)+1/2)*tau);

end do:

for i from 1 to Nx+1 do

U12[i, 1]:=fi3((i-1)*hx,((k-1)+1/2)*tau);

U12[i, Ny+1]:=fi4((i-1)*hx, ((k-1)+1/2)*tau);

end do:

# Закончили заполнять новые граничные условия

#Шаг 2 решаем ряды у, поднимаясь по х

====================================================

A_const:=b/(hy^2);

B_const:=-2*b/(hy^2)-2/tau;

C_const:=b/(hy^2);

for i from 2 to Nx do

SST:=Matrix(Ny+1, Ny+2,0);

for j from 2 to Ny do

SST[j,j-1]:=A_const;

SST[j,j]:=B_const;

SST[j,j+1]:=C_const;

SST[j, Ny+2]:=-a*(U12[i+1, j]-2*U12[i, j]+U12[i-1, j])/(hx^2) — f((i-

1)*hx,(j-1)*hy,((k-1)+1/2)*tau) -2*U12[i, j]/tau;

end do:

SST[1,1]:=C_const;

SST[Ny+1, Ny+1]:=A_const;

SST[1, Ny+2]:=U12[i,1];

SST[Ny+1, Ny+2]:=U12[i,Ny+1];

xsol:=LinearSolve(SST):

for j from 2 to Ny do

U[k+1][i,j]:=xsol[j]:

end do:

# Закончили решаеть вдоль x

# Заполняем новые граничные условия

U[k+1][i, 1]:=fi3((i-1)*hx,((k-1)+1)*tau);

U[k+1][i, Ny+1]:=fi4((i-1)*hx, ((k-1)+1)*tau);

end do:

for j from 1 to Ny+1 do

U[k+1][1, j]:=fi1((j-1)*hx,((k-1)+1)*tau);

U[k+1][Nx+1, j]:=fi2((j-1)*hy, ((k-1)+1)*tau);

end do:

# Закончили заполнять новые граничные условия

end do:

with(plots):SOLPoints:=[seq([seq([(i-1)*hx, (j-1)*hy,

U[21][i,j]],i=1..Nx+1)],j=1..Ny+1)]:

surfdata(SOLPoints, axes=normal, labels=[x,y,T]);

U_anim:=(x,y,t)-U[t][x, y]:

animate3d(U_real(x,y,t), x=1..Nx+1, y=1..Ny+1, t=1..Nt+1);

f_real:=(x,y,t)-x*y*cos(t):

plot3d(f_real(x,y,(21-1)*tau),x=0..2*Pi, y=0..2*Pi, axes=normal);

for k from 1 to Nt+1 do

for i from 1 to Nx+1 do

for j from 1 to Ny+1 do

U[k][i, j]:=f_real((i-1)*hx,(j-1)*hy,(k-1)*tau)-U[k][i, j]:

end do:

end do:

end do:

with(plots):SOLPoints1:=[seq([seq([(i-1)*hx, (j-1)*hy,

U[21][i,j]],i=1..Nx+1)],j=1..Ny+1)]:

surfdata(SOLPoints1, axes=normal, labels=[x,y,T]);

Файл «Дробных шагов.mws»

#u_t = a*u_xx + b*u_yy + f(x,y,t),

#u(0,y,t) = fi_1(y,t)=0,

#u(lx,y,t) = fi_2(y,t)=y*cos(mu*t),

#u(x,0,t) = fi_3(x,t)=0,

#u(x,ly,t) = fi_4(x,t)=x*cos(mu*t),

#u(x,y,0) = ksi(x,y)=x*y,

#f(x,y,t) = -x*y*sin(mu*t)).

restart:

#Коэффициенты теплопроводности#

a:=1:

b:=1:

mu:=1:

Nx:=20:

Ny:=20:

Nt:=30:

lx:=2*Pi:

ly:=2*Pi:

tim:=Pi/2:

fi1:=(y,t)-0.0001:

fi2:=(y,t)-y*cos(mu*t):

fi3:=(x,t)-0.0001:

fi4:=(x,t)-x*cos(mu*t):

ksi:=(x,y)-x*y:

f:=(x,y,t)—x*y*sin(mu*t):

hx:=lx/Nx:

hy:=ly/Ny:

tau:=tim/Nt:

with (LinearAlgebra):

#Задаем кубическую матрицу решения#

U:=vector(Nt+1):

for k from 1 to Nt+1 do

U[k]:=Matrix(Nx+1,Ny+1,0);

end do:

# Заполняем граничные условия

for k from 2 to Nt+1 do

for i from 1 to Nx+1 do

U[k][i, 1]:=fi3((i-1)*hx, (k-1)*tau);

U[k][i, Ny+1]:=fi4((i-1)*hx, (k-1)*tau);

end do:

for j from 1 to Ny+1 do

U[k][1, j]:=fi1((j-1)*hy, (k-1)*tau);

U[k][Nx+1, j]:=fi2((j-1)*hy, (k-1)*tau);

end do:

end do:

# Закончили заполнять граничные условия

# Заполняем начальные условия

for i from 1 to Nx+1 do

for j from 1 to Ny+1 do

U[1][i, j]:=ksi((i-1)*hx, (j-1)*hy);

end do:

end do:

# Закончили заполнять начальные условия

# Начинаем решать-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|-|

U12:=Matrix(Nx+1,Ny+1):

for k from 1 to Nt do

# Шаг 1 решаем ряды х, поднимаясь по у

====================================================

A_const:=a/(hx^2);

B_const:=-2*a/(hx^2)-1/tau;

C_const:=a/(hx^2);

for j from 2 to Ny do

SST:=Matrix(Nx+1, Nx+2,0);

for i from 2 to Nx do

SST[i,i-1]:=A_const;

SST[i,i]:=B_const;

SST[i,i+1]:=C_const;

SST[i, Nx+2]:=-f((i-1)*hx, (j-1)*hy, (k-1)*tau)/2 -U[k][i, j]/tau:

end do:

SST[1,1]:=C_const;

SST[Nx+1, Nx+1]:=A_const;

SST[1, Nx+2]:=U[k][1,j];

SST[Nx+1, Nx+2]:=U[k][Nx+1,j];

xsol:=LinearSolve(SST):

for i from 2 to Nx do

U12[i,j]:=xsol[i]:

end do:

# Закончили решаеть вдоль у

# Заполняем новые граничные условия

U12[1, j]:=fi1((j-1)*hy,((k-1)+1/2)*tau);

U12[Nx+1, j]:=fi2((j-1)*hy, ((k-1)+1/2)*tau);

end do:

for i from 1 to Nx+1 do

U12[i, 1]:=fi3((i-1)*hx,((k-1)+1/2)*tau);

U12[i, Ny+1]:=fi4((i-1)*hx, ((k-1)+1/2)*tau);

end do:

# Закончили заполнять новые граничные условия

#Шаг 2 решаем ряды у, поднимаясь по х

====================================================

A_const:=b/(hy^2):

B_const:=-2*b/(hy^2)-1/tau:

C_const:=b/(hy^2):

for i from 2 to Nx do

SST:=Matrix(Ny+1, Ny+2,0);

for j from 2 to Ny do

SST[j,j-1]:=A_const;

SST[j,j]:=B_const;

SST[j,j+1]:=C_const;

SST[j, Ny+2]:=-f((i-1)*hx, (j-1)*hy, ((k-1)+1)*tau)/2 -U12[i, j]/tau:

end do:

SST[1,1]:=C_const;

SST[Ny+1, Ny+1]:=A_const;

SST[1, Ny+2]:=U12[i,1];

SST[Ny+1, Ny+2]:=U12[i,Ny+1];

xsol:=LinearSolve(SST):

for j from 2 to Ny do

U[k+1][i,j]:=xsol[j]:

end do:

# Закончили решаеть вдоль x

# Заполняем новые граничные условия

U[k+1][i, 1]:=fi3((i-1)*hx,((k-1)+1)*tau);

U[k+1][i, Ny+1]:=fi4((i-1)*hx, ((k-1)+1)*tau);

end do:

for j from 1 to Ny+1 do

U[k+1][1, j]:=fi1((j-1)*hx,((k-1)+1)*tau);

U[k+1][Nx+1, j]:=fi2((j-1)*hy, ((k-1)+1)*tau);

end do:

# Закончили заполнять новые граничные условия

end do:

with(plots):SOLPoints:=[seq([seq([(i-1)*hx, (j-1)*hy,

U[21][i,j]],i=1..Nx+1)],j=1..Ny+1)]:

surfdata(SOLPoints, axes=normal, labels=[x,y,T]);

U_real:=(x,y,t)-U[t][x, y]:

#animate3d(U_real(x,y,t), x=1..Nx+1, y=1..Ny+1, t=1..Nt+1);

f_real:=(x,y,t)-x*y*cos(t):

plot3d(f_real(x,y,(21-1)*tau),x=0..2*Pi, y=0..2*Pi, axes=normal);

for k from 1 to Nt+1 do

for i from 1 to Nx+1 do

for j from 1 to Ny+1 do

U[k][i, j]:=f_real((i-1)*hx,(j-1)*hy,(k-1)*tau)-U[k][i, j]:

end do:

end do:

end do:

with(plots):SOLPoints1:=[seq([seq([(i-1)*hx, (j-1)*hy, U[21][i,j]],i=1..Nx+1)],j=1..Ny+1)]:

surfdata(SOLPoints1, axes=normal, labels=[x,y,T]);

Файл « diff.mws»

#u_t = a*u_xx + b*u_yy + f(x,y,t),

#u(0,y,t) = fi_1(y,t)=0,

#u(lx,y,t) = fi_2(y,t)=y*cos(mu*t),

#u(x,0,t) = fi_3(x,t)=0,

#u(x,ly,t) = fi_4(x,t)=x*cos(mu*t),

#u(x,y,0) = ksi(x,y)=0,

#f(x,y,t) = x*y*sin(mu*t)).

restart:

#Коэффициенты теплопроводности#

a:=1:

b:=1:

mu:=1:

Nx:=20:

Ny:=20:

Nt:=30:

lx:=2*Pi:

ly:=2*Pi:

tim:=Pi/2:

fi1:=(y,t)-0.0001:

fi2:=(y,t)-y*cos(mu*t):

fi3:=(x,t)-0.0001:

fi4:=(x,t)-x*cos(mu*t):

ksi:=(x,y)-x*y:

f:=(x,y,t)—x*y*sin(mu*t):

hx:=lx/Nx:

hy:=ly/Ny:

tau:=tim/Nt:

with (LinearAlgebra):

#Задаем кубическую матрицу решения#

U:=vector(Nt+1):

for k from 1 to Nt+1 do

U[k]:=Matrix(Nx+1,Ny+1,0);

end do:

# Заполняем граничные условия

for k from 2 to Nt+1 do

for i from 1 to Nx+1 do

U[k][i, 1]:=fi3((i-1)*hx, (k-1)*tau);

U[k][i, Ny+1]:=fi4((i-1)*hx, (k-1)*tau);

end do:

for j from 1 to Ny+1 do

U[k][1, j]:=fi1((j-1)*hy, (k-1)*tau);

U[k][Nx+1, j]:=fi2((j-1)*hy, (k-1)*tau);

end do:

end do:

# Закончили заполнять граничные условия

# Заполняем начальные условия

for i from 1 to Nx+1 do

for j from 1 to Ny+1 do

U[1][i, j]:=ksi((i-1)*hx, (j-1)*hy);

end do:

end do:

# Закончили заполнять начальные условия

# Начинаем решать вариант переменных направлений-|-|-||-|-|-|-|-|-|-|-|-|-|-|

U12:=Matrix(Nx+1,Ny+1):

for k from 1 to Nt do

# Шаг 1 решаем ряды х, поднимаясь по у

====================================================

A_const:=a/(hx^2);

B_const:=-2*a/(hx^2)-2/tau;

C_const:=a/(hx^2);

for j from 2 to Ny do

SST:=Matrix(Nx+1, Nx+2,0);

for i from 2 to Nx do

SST[i,i-1]:=A_const;

SST[i,i]:=B_const;

SST[i,i+1]:=C_const;

SST[i, Nx+2]:=-b*(U[k][i, j+1]-2*U[k][i, j]+U[k][i, j-1])/(hy^2) — f((i-

1)*hx,(j-1)*hy,((k-1)+1/2)*tau) -2*U[k][i, j]/tau;

end do:

SST[1,1]:=C_const;

SST[Nx+1, Nx+1]:=A_const;

SST[1, Nx+2]:=U[k][1,j];

SST[Nx+1, Nx+2]:=U[k][Nx+1,j];

xsol:=LinearSolve(SST):

for i from 2 to Nx do

U12[i,j]:=xsol[i]:

end do:

# Закончили решаеть вдоль у

# Заполняем новые граничные условия

U12[1, j]:=fi1((j-1)*hy,((k-1)+1/2)*tau);

U12[Nx+1, j]:=fi2((j-1)*hy, ((k-1)+1/2)*tau);

end do:

for i from 1 to Nx+1 do

U12[i, 1]:=fi3((i-1)*hx,((k-1)+1/2)*tau);

U12[i, Ny+1]:=fi4((i-1)*hx, ((k-1)+1/2)*tau);

end do:

# Закончили заполнять новые граничные условия

#Шаг 2 решаем ряды у, поднимаясь по х

====================================================

A_const:=b/(hy^2);

B_const:=-2*b/(hy^2)-2/tau;

C_const:=b/(hy^2);

for i from 2 to Nx do

SST:=Matrix(Ny+1, Ny+2,0);

for j from 2 to Ny do

SST[j,j-1]:=A_const;

SST[j,j]:=B_const;

SST[j,j+1]:=C_const;

SST[j, Ny+2]:=-a*(U12[i+1, j]-2*U12[i, j]+U12[i-1, j])/(hx^2) — f((i-

1)*hx,(j-1)*hy,((k-1)+1/2)*tau) -2*U12[i, j]/tau;

end do:

SST[1,1]:=C_const;

SST[Ny+1, Ny+1]:=A_const;

SST[1, Ny+2]:=U12[i,1];

SST[Ny+1, Ny+2]:=U12[i,Ny+1];

xsol:=LinearSolve(SST):

for j from 2 to Ny do

U[k+1][i,j]:=xsol[j]:

end do:

# Закончили решаеть вдоль x

# Заполняем новые граничные условия

U[k+1][i, 1]:=fi3((i-1)*hx,((k-1)+1)*tau);

U[k+1][i, Ny+1]:=fi4((i-1)*hx, ((k-1)+1)*tau);

end do:

for j from 1 to Ny+1 do

U[k+1][1, j]:=fi1((j-1)*hx,((k-1)+1)*tau);

U[k+1][Nx+1, j]:=fi2((j-1)*hy, ((k-1)+1)*tau);

end do:

# Закончили заполнять новые граничные условия

end do:

U1:=vector(Nt+1):

for k from 1 to Nt+1 do

U1[k]:=Matrix(Nx+1,Ny+1,0);

end do:

for k from 1 to Nt+1 do

for i from 1 to Nx+1 do

for j from 1 to Ny+1 do

U1[k][i, j]:=U[k][i, j]:

end do:

end do:

end do:

#u_t = a*u_xx + b*u_yy + f(x,y,t),

#u(0,y,t) = fi_1(y,t)=0,

#u(lx,y,t) = fi_2(y,t)=y*cos(mu*t),

#u(x,0,t) = fi_3(x,t)=0,

#u(x,ly,t) = fi_4(x,t)=x*cos(mu*t),

#u(x,y,0) = ksi(x,y)=x*y,

#f(x,y,t) = -x*y*sin(mu*t)).

# Начинаем решать вариант дробных шагов-|-|-|-|-|-|-|-|-|-|-|-|-|—|-|-|-|-|-|-|-|-|

U12:=Matrix(Nx+1,Ny+1):

for k from 1 to Nt do

# Шаг 1 решаем ряды х, поднимаясь по у

====================================================

A_const:=a/(hx^2);

B_const:=-2*a/(hx^2)-1/tau;

C_const:=a/(hx^2);

for j from 2 to Ny do

SST:=Matrix(Nx+1, Nx+2,0);

for i from 2 to Nx do

SST[i,i-1]:=A_const;

SST[i,i]:=B_const;

SST[i,i+1]:=C_const;

SST[i, Nx+2]:=-f((i-1)*hx, (j-1)*hy, (k-1)*tau)/2 -U[k][i, j]/tau:

end do:

SST[1,1]:=C_const;

SST[Nx+1, Nx+1]:=A_const;

SST[1, Nx+2]:=U[k][1,j];

SST[Nx+1, Nx+2]:=U[k][Nx+1,j];

xsol:=LinearSolve(SST):

for i from 2 to Nx do

U12[i,j]:=xsol[i]:

end do:

# Закончили решаеть вдоль у

# Заполняем новые граничные условия

U12[1, j]:=fi1((j-1)*hy,((k-1)+1/2)*tau);

U12[Nx+1, j]:=fi2((j-1)*hy, ((k-1)+1/2)*tau);

end do:

for i from 1 to Nx+1 do

U12[i, 1]:=fi3((i-1)*hx,((k-1)+1/2)*tau);

U12[i, Ny+1]:=fi4((i-1)*hx, ((k-1)+1/2)*tau);

end do:

# Закончили заполнять новые граничные условия

#Шаг 2 решаем ряды у, поднимаясь по х

====================================================

A_const:=b/(hy^2):

B_const:=-2*b/(hy^2)-1/tau:

C_const:=b/(hy^2):

for i from 2 to Nx do

SST:=Matrix(Ny+1, Ny+2,0);

for j from 2 to Ny do

SST[j,j-1]:=A_const;

SST[j,j]:=B_const;

SST[j,j+1]:=C_const;

SST[j, Ny+2]:=-f((i-1)*hx, (j-1)*hy, ((k-1)+1)*tau)/2 -U12[i, j]/tau:

end do:

SST[1,1]:=C_const;

SST[Ny+1, Ny+1]:=A_const;

SST[1, Ny+2]:=U12[i,1];

SST[Ny+1, Ny+2]:=U12[i,Ny+1];

xsol:=LinearSolve(SST):

for j from 2 to Ny do

U[k+1][i,j]:=xsol[j]:

end do:

# Закончили решать вдоль x

# Заполняем новые граничные условия

U[k+1][i, 1]:=fi3((i-1)*hx,((k-1)+1)*tau);

U[k+1][i, Ny+1]:=fi4((i-1)*hx, ((k-1)+1)*tau);

end do:

for j from 1 to Ny+1 do

U[k+1][1, j]:=fi1((j-1)*hx,((k-1)+1)*tau);

U[k+1][Nx+1, j]:=fi2((j-1)*hy, ((k-1)+1)*tau);

end do:

# Закончили заполнять новые граничные условия

end do:

for k from 1 to Nt+1 do

for i from 1 to Nx+1 do

for j from 1 to Ny+1 do

U[k][i, j]:=U1[k][i, j]-U[k][i, j]:

end do:

end do:

end do:

with(plots):SOLPoints:=[seq([seq([(i-1)*hx, (j-1)*hy,

U[21][i,j]],i=1..Nx+1)],j=1..Ny+1)]:

surfdata(SOLPoints, axes=normal, labels=[x,y,T]);