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

% 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;

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

%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);

%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);