Ирина Эланс

Автор который поможет с любыми образовательными и учебными заданиями

Численное решение уравнений теплопроводности методом разностных схем

ФЕДЕРАЛЬНОЕ АГЕНТСТВО ПО ОБРАЗОВАНИЮ

Филиал федерального государственного автономного образовательного учреждения высшего профессионального образования «Казанский (Приволжский) федеральный университет» в г. Набережные Челны.

ФАКУЛЬТЕТ ПРИКЛАДНОЙ МАТЕМАТИКИ
И ИНФОРМАЦИОННЫХ ТЕХНОЛОГИЙ

КАФЕДРА ПРИКЛАДНОЙ МАТЕМАТИКИ И ИНФОРМАТИКИ

Специальность: 010501.65 – Прикладная математика и информатика

КУРСОВАЯ РАБОТА

IX СЕМЕСТР

ТЕМА: «Решение нелинейных уравнений переноса тепла методом разностных схем»

Дисциплина: численные методы

Выполнил

студент Токмаков А.М.

группа 4606 курс 5

Научный руководитель

Марданшин Р. Г.

к. ф.–м. н., доцент

Члены комиссии по защите курсовой работы

_________________________________________

^{Ф.И.О.
подпись}

_________________________________________

^{Ф.И.О.
подпись}

_________________________________________

^{Ф.И.О.
подпись}

Оценка________________

Дата ________________

Набережные Челны

2010

Оглавление

Аннотация

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

Введение

Численное решение прикладных задач всегда интересовало математиков. Крупнейшие представители прошлого сочетали в своих исследованиях изучение явлений природы, получение их математического описания, как иногда говорят, математической модели явления, и его исследование. Анализ усложненных моделей потребовал создание специальных, как правило, численных или асимптотических методов решения задач. Названия некоторых из таких методов – методы Ньютона, Эйлера, Лобачевского, Гаусса, Чебышева, Эрмита, Крылова – свидетельствуют о том, что их разработкой занимались крупнейшие ученые своего времени.

Настоящее время характерно резким расширением приложений математики, во многим связанным с созданием и развитием средств вычислительной техники. В результате появления ЭВМ (электронно-вычислительных машин, или как часто говорят, компьютеров) с программным управлением менее чем за 50 лет скорость выполнения арифметических операций возросла от 0.1 операции в секунду при ручном расчете до 10¹² операций на современных серийных ЭВМ, т.е. примерно в 10¹³ раз.

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

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

Цель работы: разработать разностные схемы, позволяющие решать задачу нелинейные уравнения переноса тепла. В качестве среды разработки был выбран пакет matlab 6.5.

Математическое описание основных понятий и разностных схем

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

Рассмотрим сначала некоторые необходимые понятия из теории сеток:

Пусть имеется пространство , где - функция непрерывного аргумента . На отрезке введем конечное множество точек , которое назовем сеткой. Точки , будем называть узлами сетки . Множество без узлов и будем обозначать . Если расстояние между соседними узлами постоянно (не зависит от i), для всех , то сетку называют равномерной (с шагом h), в противном случае – неравномерной. Вместо функции , определенной для всех , будем рассматривать сеточную функцию , целочисленного аргумента или узла сетки , а заменим конечномерным (размерностью N+1) пространством сеточных функций. Очевидно, что сеточную функцию можно рассматривать как вектор [1].

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

Рассмотрим первую краевую задачу для нелинейного уравнения теплопроводности в следующем виде [3]:

(1)

, ,

, .

в области . Требуется найти непрерывное в решение .

Разностные схемы [3]

В области введем сетку

Заменяем производную по t разностными выражением [2]

с шагом по t. Дифференциальное выражение при каждом фиксированном t аппроксимируем следующим образом [4]. Сначала берем левостороннюю разностную аппроксимацию производной по x, а после правостороннюю.

(2)

В граничных условиях заменяем непрерывную функцию дискретной:

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

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

(3)

Эта схема абсолютно устойчива и имеет первый порядок аппроксимации по и второй – по , решение находится методом прогонки. Перепишем (3) в другом виде

(4)

Для реализации этой схемы необходимо применить тот или иной итерационный метод. Например такой:

(5)

Здесь s – номер итерации. Как видим, нелинейные коэффициенты берутся с предыдущей итерации, а качестве начального приближения для выбирается . Это начальное приближение тем лучше, чем меньше шаг . Число итераций M задается из соображения точности. В задачах с гладкими коэффициентами при часто бывает достаточно провести две - три итерации. Значения на новой итерации находятся методом прогонки. При M=1 итерационный метод (5) совпадает с разностной схемой (4).

Для приближенного решения нелинейной схемы (4) применяются также схемы предиктор-корректор второго порядка точности , аналогичные методу Рунге-Кутта для обыкновенных дифференциальных уравнений. Здесь переход со слоя k на k+1 осуществляется в два этапа. На первом этапе решается линейная неявная система уравнений

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

Библиотека функций, разработки

В работе разработаны следующие функции для решения неоднородного одномерного уравнения теплопроводности:

TeploProvodNotLine1(X,N,T,K,k_U,f_U,U0,U1,U2) – находит приближенное решение поставленной задачи при помощи линейной неявной схемы.
TeploProvodNotLine(X,N,T,K,k_U,f_U,U0,U1,U2) – находит приближенное решение поставленной задачи при помощи схемы предиктор-корректор.
TeploProvodNotLine2(X,N,T,K,k_U,f_U,U0,U1,U2,E) – находит приближенное решение поставленной задачи при помощи нелинейной неявной схемы

где X – длина стержня, N – количество узлов сетки, по пространственной переменной, T - время, K – число узлов по времени, k_U – функция нелинейности, f_U – функция источника, U0 – начальное распределение температуры, U1 – первое граничное условие, U2 – второе граничное условие, E – точность вычислений.

Функции k_U. f_U, U0, U1, U2 задаются в виде указателей на m-файлы. Для обеспечения универсальности.

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

Методы объединяет экранная форма, написанная в среде matlab 6.5 [6]:

На ней присутствуют следующие элементы:

Три радио-кнопки, позволяющие выбрать один из разработанных трех методов.
Параметры задачи (длинна стержня, время, количество узлов сетки по пространственной координате, количество узлов по времени, k_U -функция нелинейности, f_U – функция источника, U0 – начальное распределение температуры в стержне, U1,U2 – граничные условия первого рода, точность (используется для нелинейной неявной схемы)).
Кнопка «Считать», обработка события нажатия которой осуществляет вычисления и выводит график поверхности распределения температуры.

Тестирование

Протестируем разработанные методы на примере уравнения Фишера [5].

. Примем следующие параметры разработанных функций:

X=1 длина стержня;

N=25 – кол-во узлов сетки по пространственной переменной;

T=0.1 - время;

K=10 кол-во узлов сетки по временной переменной;

f(U)=0.45*U*(1-U)– функция граничных условий;

k(U)=1;

E =0.001 – точность;

U0= sin(X+0.45) – начальное распределение температуры;

U1= cos(T+0.45) – первое граничное условие;

U2= 6*T+0.9 – второе граничное условие;

Выберем линейную неявную схему

Рисунок 1

Результат

Рисунок 2

Теперь с теми же параметрами произведем расчеты с помощью нелинейной неявной схемы

Рисунок 3

Полученный результат

Рисунок 4

Наша функция , таким образом в нашем случае мы имеем дело с гладкими коэффициентами. И здесь достаточно 2-3 итераций для получения решения. Как видно график мало чем отличается от линейной явной схемы.

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

Результат линейной явной схемы

Рисунок 5

Для нелинейной схемы зададим точность повыше E=0.0000001.

Результат нелинейной неявной схемы

Рисунок 6

Протестируем схему предиктор-корректор при исходных условиях

Рисунок 7

Результат имеет некоторое отличие от предыдущих двух

Рисунок 8

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

Проверим поведения схемы при негладких коэффициентах

Рисунок 9

Результат имеет серьезные отличия от предыдущих схем.

Вывод

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

Заключение

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

Также были усвоены и закреплены навыки:

Создание графических экранных форм в среде matlab 6.5.
Программирования в среде matlab.
Использование указателей на функции.
Разработки численных методов для нелинейных уравнений переноса тепла.

Список использованной литературы

Самарский А.А. Введение в численные методы. Учебное пособие для вузов. 3-е изд., стер. –СПб.: Издательство «Лань», 2005. – 288 с.: ил. – (Учебники для вузов. Специальная литература).
Бахвалов Н.С., Жидков Н. П., Кобельков Г. М. Численные методы – 3-е изд, доп., и перераб. – М.: БИНОМ. Лаборатория знаний, 2004. – 636 с., илл.
Самарский А.А., Гулин А.В., Численные методы: Учеб. пособие для ВУЗов. – М.: Наука. Гл. ред. Физ-мат. лит., 1989 – 432 с.
Турчак Л. И., Плотников П. В. Основы численных методов: Учеб. пособие. – 2-е изд., перераб. И доп. – М.: ФИЗМАТЛИТ, 2005. -304 с.
http://eqworld.ipmnet.ru/ru/solutions/npde/npde-toc1.htm - электронная книга
http://matlab.exponenta.ru/gui/index.php - электронная книга.

Приложения

GUI.m

function varargout = GUI(varargin)

% GUI M-file for GUI.fig

% GUI, by itself, creates a new GUI or raises the existing

% singleton*.

% H = GUI returns the handle to a new GUI or the handle to

% the existing singleton*.

% GUI('CALLBACK',hObject,eventData,handles,...) calls the local

% function named CALLBACK in GUI.M with the given input arguments.

% GUI('Property','Value',...) creates a new GUI or raises the

% existing singleton*. Starting from the left, property value pairs are

% applied to the GUI before GUI_OpeningFunction gets called. An

% unrecognized property name or invalid value makes property application

% stop. All inputs are passed to GUI_OpeningFcn via varargin.

% *See GUI Options on GUIDE's Tools menu. Choose "GUI allows only one

% instance to run (singleton)".

% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help GUI

% Last Modified by GUIDE v2.5 04-Nov-2010 23:32:14

% Begin initialization code - DO NOT EDIT

gui_Singleton = 1;

gui_State = struct('gui_Name', mfilename, ...

'gui_Singleton', gui_Singleton, ...

'gui_OpeningFcn', @GUI_OpeningFcn, ...

'gui_OutputFcn', @GUI_OutputFcn, ...

'gui_LayoutFcn', [] , ...

'gui_Callback', []);

if nargin & isstr(varargin{1})

gui_State.gui_Callback = str2func(varargin{1});

end

if nargout

[varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});

else

gui_mainfcn(gui_State, varargin{:});

end

% End initialization code - DO NOT EDIT

% --- Executes just before GUI is made visible.

function GUI_OpeningFcn(hObject, eventdata, handles, varargin)

% This function has no output args, see OutputFcn.

% hObject handle to figure

% eventdata reserved - to be defined in a future version of MATLAB

% handles structure with handles and user data (see GUIDATA)

% varargin command line arguments to GUI (see VARARGIN)

% Choose default command line output for GUI

handles.output = hObject;

% Update handles structure

guidata(hObject, handles);

set(handles.radiobutton4,'Value',1);

set(handles.radiobutton5,'Value',0);

set(handles.radiobutton6,'Value',0);

set(handles.edit11,'Enable','inactive');

% UIWAIT makes GUI wait for user response (see UIRESUME)

% uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.

function varargout = GUI_OutputFcn(hObject, eventdata, handles)

% varargout cell array for returning output args (see VARARGOUT);

% hObject handle to figure

% eventdata reserved - to be defined in a future version of MATLAB

% handles structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure

varargout{1} = handles.output;

% --- Executes on button press in pushbutton1.

function pushbutton1_Callback(hObject, eventdata, handles)

X=str2num(get(handles.edit2,'String'));

N=str2num(get(handles.edit3,'String'));

T=str2num(get(handles.edit4,'String'));

K=str2num(get(handles.edit5,'String'));

f_U=str2num(get(handles.edit6,'String'));

k_U=str2num(get(handles.edit7,'String'));

U0=str2num(get(handles.edit8,'String'));

U1=str2num(get(handles.edit9,'String'));

U2=str2num(get(handles.edit10,'String'));

key1=get(handles.radiobutton4,'Value');

key2=get(handles.radiobutton5,'Value');

if (key1==1) %Если выбрана линейная неявная схема

TeploprovodNotLine1(X,N,T,K,k_U,f_U,U0,U1,U2);

else

if(key2==1) %Если выбрана схема предиктор-корректор

TeploprovodNotLine(X,N,T,K,k_U,f_U,U0,U1,U2);

else

E=str2num(get(handles.edit11,'String'));

%Выбрана нелинейная неявная схема

TeploprovodNotLine2(X,N,T,K,k_U,f_U,U0,U1,U2,E);

end;

% hObject handle to pushbutton1 (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB

% handles structure with handles and user data (see GUIDATA)

% --- Executes on button press in radiobutton4.

function radiobutton4_Callback(hObject, eventdata, handles)

set(handles.radiobutton4,'Value',1);

set(handles.radiobutton5,'Value',0);

set(handles.radiobutton6,'Value',0);

set(handles.edit11,'Enable','inactive');

% hObject handle to radiobutton4 (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB

% handles structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of radiobutton4

% --- Executes on button press in radiobutton5.

function radiobutton5_Callback(hObject, eventdata, handles)

set(handles.radiobutton4,'Value',0);

set(handles.radiobutton5,'Value',1);

set(handles.radiobutton6,'Value',0);

set(handles.edit11,'Enable','inactive');

% hObject handle to radiobutton5 (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB

% handles structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of radiobutton5

% --- Executes during object creation, after setting all properties.

function edit1_CreateFcn(hObject, eventdata, handles)

% hObject handle to edit1 (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB

% handles empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.

% See ISPC and COMPUTER.

if ispc

set(hObject,'BackgroundColor','white');

else

set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));

end

function edit1_Callback(hObject, eventdata, handles)

% hObject handle to edit1 (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB

% handles structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit1 as text

% str2double(get(hObject,'String')) returns contents of edit1 as a double

% --- Executes on button press in radiobutton6.

function radiobutton6_Callback(hObject, eventdata, handles)

set(handles.radiobutton4,'Value',0);

set(handles.radiobutton5,'Value',0);

set(handles.radiobutton6,'Value',1);

set(handles.edit11,'Enable','on');

% hObject handle to radiobutton6 (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB

% handles structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of radiobutton6

% --- Executes during object creation, after setting all properties.

function edit2_CreateFcn(hObject, eventdata, handles)

% hObject handle to edit2 (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB

% handles empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.

% See ISPC and COMPUTER.

if ispc

set(hObject,'BackgroundColor','white');

else

set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));

end

function edit2_Callback(hObject, eventdata, handles)

% hObject handle to edit2 (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB

% handles structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit2 as text

% str2double(get(hObject,'String')) returns contents of edit2 as a double

% --- Executes during object creation, after setting all properties.

function edit3_CreateFcn(hObject, eventdata, handles)

% hObject handle to edit3 (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB

% handles empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.

% See ISPC and COMPUTER.

if ispc

set(hObject,'BackgroundColor','white');

else

set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));

end

function edit3_Callback(hObject, eventdata, handles)

% hObject handle to edit3 (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB

% handles structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit3 as text

% str2double(get(hObject,'String')) returns contents of edit3 as a double

% --- Executes during object creation, after setting all properties.

function edit4_CreateFcn(hObject, eventdata, handles)

% hObject handle to edit4 (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB

% handles empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.

% See ISPC and COMPUTER.

if ispc

set(hObject,'BackgroundColor','white');

else

set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));

end

function edit4_Callback(hObject, eventdata, handles)

% hObject handle to edit4 (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB

% handles structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit4 as text

% str2double(get(hObject,'String')) returns contents of edit4 as a double

% --- Executes during object creation, after setting all properties.

function edit5_CreateFcn(hObject, eventdata, handles)

% hObject handle to edit5 (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB

% handles empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.

% See ISPC and COMPUTER.

if ispc

set(hObject,'BackgroundColor','white');

else

set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));

end

function edit5_Callback(hObject, eventdata, handles)

% hObject handle to edit5 (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB

% handles structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit5 as text

% str2double(get(hObject,'String')) returns contents of edit5 as a double

% --- Executes during object creation, after setting all properties.

function edit6_CreateFcn(hObject, eventdata, handles)

% hObject handle to edit6 (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB

% handles empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.

% See ISPC and COMPUTER.

if ispc

set(hObject,'BackgroundColor','white');

else

set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));

end

function edit6_Callback(hObject, eventdata, handles)

% hObject handle to edit6 (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB

% handles structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit6 as text

% str2double(get(hObject,'String')) returns contents of edit6 as a double

% --- Executes during object creation, after setting all properties.

function edit7_CreateFcn(hObject, eventdata, handles)

% hObject handle to edit7 (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB

% handles empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.

% See ISPC and COMPUTER.

if ispc

set(hObject,'BackgroundColor','white');

else

set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));

end

function edit7_Callback(hObject, eventdata, handles)

% hObject handle to edit7 (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB

% handles structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit7 as text

% str2double(get(hObject,'String')) returns contents of edit7 as a double

% --- Executes during object creation, after setting all properties.

function edit8_CreateFcn(hObject, eventdata, handles)

% hObject handle to edit8 (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB

% handles empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.

% See ISPC and COMPUTER.

if ispc

set(hObject,'BackgroundColor','white');

else

set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));

end

function edit8_Callback(hObject, eventdata, handles)

% hObject handle to edit8 (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB

% handles structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit8 as text

% str2double(get(hObject,'String')) returns contents of edit8 as a double

% --- Executes during object creation, after setting all properties.

function edit9_CreateFcn(hObject, eventdata, handles)

% hObject handle to edit9 (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB

% handles empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.

% See ISPC and COMPUTER.

if ispc

set(hObject,'BackgroundColor','white');

else

set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));

end

function edit9_Callback(hObject, eventdata, handles)

% hObject handle to edit9 (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB

% handles structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit9 as text

% str2double(get(hObject,'String')) returns contents of edit9 as a double

% --- Executes during object creation, after setting all properties.

function edit10_CreateFcn(hObject, eventdata, handles)

% hObject handle to edit10 (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB

% handles empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.

% See ISPC and COMPUTER.

if ispc

set(hObject,'BackgroundColor','white');

else

set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));

end

function edit10_Callback(hObject, eventdata, handles)

% hObject handle to edit10 (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB

% handles structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit10 as text

% str2double(get(hObject,'String')) returns contents of edit10 as a double

% --- Executes during object creation, after setting all properties.

function edit11_CreateFcn(hObject, eventdata, handles)

% hObject handle to edit11 (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB

% handles empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.

% See ISPC and COMPUTER.

if ispc

set(hObject,'BackgroundColor','white');

else

set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));

end

function edit11_Callback(hObject, eventdata, handles)

% hObject handle to edit11 (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB

% handles structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit11 as text

% str2double(get(hObject,'String')) returns contents of edit11 as a double

TeploprovodNotLine.m

function [x,t,Ymain]=TeploProvodNotLine(l, N, T, K, k_U, f_U, U0, U1, U2) %Схема предиктор-корректор

%l - длинна стержня

%N - Разбиение сетки по пространственной координате

%T - Разбиение сетки по временной координате

%k_U - функция k(U)

%f_U - функция f(U)

%U0 - Начальное распределение температуры

%U1 - Функция подачи тепла на начало стержня

%U2 - Функция подачи тепла на конец стержня

%-----------------------------------------

%создание сетки

h=l/N;

tau=T/K;

x=[0:h:h*N];

t=[0:tau:tau*K];

coef=0.5*tau/h^2;

%-----------------------------------------

Y05=zeros(1,N+1); % предикторы

Ymain=zeros(K+1,N+1); % корректоры - матрица полученных значений.

Ymain(1,:)=feval(U0,x); % Вставка значений начальных условий в первую строку, где U(x,0)=f(x);

Ymain(:,1)=feval(U1,t'); % Вставка значений начальных условий в первый столбец, где U(0,t)=f1(t);

Ymain(:,N+1)=feval(U2,t'); % Вставка значений начальных условий в последний столбец, где U(1,t)=f2(t);

U=zeros(N-1,N-1);

UpDiag=zeros(1,N);

a=zeros(1,N);

B=zeros(1,N-1);

PredSolution=zeros(1, N+1);

for i=1:K % цикл по слоям времени

%Вычисление предикторов

PredSolution(1,1)=feval(U1,t(i)+0.5*tau);

PredSolution(1,N+1)=feval(U2,t(i)+0.5*tau);

%a=feval(k_U,Ymain(i,2:N+1)); % функция a(y) по текущему слою для предикторов {Переопределить!!!}

for j=1:N

%display(Ymain(i,j+1));

%lp=feval(k_U,Ymain(i,j+1));

%display(Ymain(i,j));

%rp=feval(k_U,Ymain(i,j));

a(j)=0.5*(feval(k_U,Ymain(i,j+1))+feval(k_U,Ymain(i,j))); %Всего N элементов, крайние подут в правые части

end;

%составление основной матрицы для С.Л.А.У. предикторов

UpDiag=(0.5*tau/h^2).*a;

U=diag(UpDiag(1,2:N-1),1)+diag(UpDiag(1,2:N-1),-1);

for j=1:N-1

U(j,j)=-(1+((0.5*tau)/h^2)*(a(j)+a(j+1)));

end

% Составление правых частей для С.Л.А.У. предикторов

B=-0.5*tau*feval(f_U,Ymain(i,2:N))-Ymain(i,2:N);

B(1,1)=B(1,1)-UpDiag(1,1)*PredSolution(1,1);

B(1,N-1)=B(1,N-1)-UpDiag(1,N)*PredSolution(1,N+1);

X=U\B';

PredSolution(1,2:N)=X';

%Вычисление корректоров

%a=feval(k_U, PredSolution'); % Вычисление функции a(y) от вычисленных предикторов

for j=1:N

a(j)=0.5*(feval(k_U,PredSolution(j+1))+feval(k_U,PredSolution(j))); %Всего N элементов, крайние подут в правые части

end;

UpDiag=(0.5*tau/h^2).*a; % Нахождение главной диагонали для матрицы корректоров

U=diag(UpDiag(1,2:N-1),1)+diag(UpDiag(1,2:N-1),-1);

for j=1:N-1

U(j,j)=-(1+((0.5*tau)/h^2)*(a(j)+a(j+1))); % нахождение главной диагонали для матрицы корректоров

end

%Составление правых частей для корректоров

for j=1:N-1

q1=-0.5*tau*a(j)*Ymain(i,j);

q2=(-1+(0.5*tau/(h^2))*(a(j+1)+a(j)))*Ymain(i,j+1);

q3=(-0.5*tau/h^2)*a(j+1)*Ymain(i,j+2);

q4=-tau*feval(k_U,PredSolution(j+1));

B(1,j)=-(0.5*tau/h^2)*a(j)*Ymain(i,j)+(-1+(0.5*tau/h^2)*(a(j+1)+a(j)))*Ymain(i,j+1)-(0.5*tau/h^2)*a(j+1)*Ymain(i,j+2)-tau*feval(f_U,PredSolution(j+1));

end

B(1,1)=B(1,1)-(0.5*tau/h^2)*a(1)*Ymain(i+1,1);

B(1,N-1)=B(1,N-1)-(0.5*tau/h^2)*a(N)*Ymain(i+1,N+1);

Solution=U\B';

Ymain(i+1,2:N)=Solution';

%for n=1:N

%end;

%for n=1:N

%end;

end;

figure;

surf(x,t,Ymain);

TeploProvodNotLine1.m

function [x,t,Ymain]=TeploprovodNotLine1(l, N, T, K, k_U, f_U, U0, U1, U2) % линейная неявная схема

%l - длинна стержня

%N - Разбиение сетки по пространственной координате

%T - Разбиение сетки по временной координате

%k_U - функция k(U)

%f_U - функция f(U)

%U0 - Начальное распределение температуры

%U1 - Функция подачи тепла на начало стержня

%U2 - Функция подачи тепла на конец стержня

%-----------------------------------------

%создание сетки

h=l/N;

tau=T/K;

x=[0:h:h*N];

t=[0:tau:tau*K];

coef=0.5*tau/h^2;

%-----------------------------------------

Ymain=zeros(K+1,N+1); % корректоры - матрица полученных значений.

Ymain(1,:)=feval(U0,x); % Вставка значений начальных условий в первую строку, где U(x,0)=f(x);

Ymain(:,1)=feval(U1,t'); % Вставка значений начальных условий в первый столбец, где U(0,t)=f1(t);

Ymain(:,N+1)=feval(U2,t'); % Вставка значений начальных условий в последний столбец, где U(1,t)=f2(t);

U=zeros(N-1,N-1);

for i=1:K

%заполняем ai

for j=1:N

a(j)=0.5*(feval(k_U,Ymain(i,j+1))+feval(k_U,Ymain(i,j))); %Всего N элементов, крайние подут в правые части

end;

%составление основной матрицы для С.Л.А.У.

%диагонали, расположенные вышеи ниже главной

UpDiag=(tau/h^2).*a;

U=diag(UpDiag(1,2:N-1),1)+diag(UpDiag(1,2:N-1),-1);

%Главная диагональ

for j=1:N-1

U(j,j)=-(1+(tau/h^2)*(a(j)+a(j+1)));

end;

% Составление правых частей для С.Л.А.У.

B=-tau*feval(f_U,Ymain(i,2:N))-Ymain(i,2:N);

B(1,1)=B(1,1)-(tau/h^2)*a(1)*Ymain(i+1,1);

B(1,N-1)=B(1,N-1)-(tau/h^2)*a(N)*Ymain(i+1,N+1);

Solution=U\B';

Ymain(i+1,2:N)=Solution';

end;

figure;

surf(x,t,Ymain);