Commit b9b11430 authored by Guilherme Tegoni Goedert's avatar Guilherme Tegoni Goedert
Browse files

added some code from the 2016 course

parent 5df0b18f
function dw = Lorentz(t,w,sigma,rho,beta)
dw(1) = sigma*(w(2) - w(1));
dw(2) = w(1)*(rho - w(3)) - w(2);
dw(3) = w(1)*w(2) - beta*w(3);
end
\ No newline at end of file
% main
% Initial conditions
t0 = 0;
w0 = [1,1,1];
% Maximum time of integration
T = 100;
% Step definitions
h0 = 0.01; % initial step size
h = h0; % used step size
hmin = (1e-5)*h0; % minimum step size allowed
% Model parameters and definitions
sigma = 10;
rho = 28;
beta = 8/3;
F = @(s,w) Lorenz(s,w,sigma,rho,beta);
% Poincare Section definition
Section = @ (w) w(3) - 15;
tsec = [];
wsec = [];
% Result allocation
t = zeros(100000,1);
w = zeros(100000,max(size(w0)));
t(1) = t0;
w(1,:) = w0;
counter = 2;
% Initiate figure object
figure;
az = 0; el = 0; % viewpoint definition
%%%%%% ITERATE
while t(counter) <= T
% iterate a step
[t(counter), w(counter,:)] = RK4(t(counter-1),w(counter-1,:),h,F);
% evaluete if solutuion passes through the poincare section
if or(sign(Section(w(counter,:))) > sign(Section(w(counter-1,:))), sign(Section(w(counter,:))) < sign(Section(w(counter-1,:))))
h = h/2;
if h <= hmin
h = h0;
tsec = [tsec;t(counter)];
wsec = [wsec;w(counter,:)];
counter = counter + 1;
end
else
% plot result
plot3(w(1:counter,1),w(1:counter,2),w(1:counter,3))
xlabel('x'); ylabel('y'); zlabel('z')
hold on;
% plot points that pass through Poincare section
if max(size(wsec))>0
plot3(wsec(:,1),wsec(:,2),wsec(:,3),'o');
end
% plot current point in trajectory
plot3(w(counter,1),w(counter,2),w(counter,3),'*')
hold off
% rotate and view solution
az = mod(az + 0.5,360);
view(az,el)
pause(0.005);
%iterate
counter = counter + 1;
end
end
% eliminate unused entries in solution vector
t(counter+1:end) = [];
w(counter+1,:) = [];
% main
t0 = 0;
w0 = [1,1,1];
T = 100;
h = 0.01;
sigma = 10;
rho = 28;
beta = 8/3;
F = @(s,w) Lorentz(s,w,sigma,rho,beta);
t = zeros(2000,1);
w = zeros(2000,max(size(w0)));
t(1) = t0;
w(1,:) = w0;
counter = 1;
figure;
az = 0; el = 45;
while t <= T
counter = counter + 1;
[t(counter), w(counter,:)] = RK4(t(counter-1),w(counter-1,:),h,F);
plot3(w(1:counter,1),w(1:counter,2),w(1:counter,3))
hold on;
plot3(w(counter,1),w(counter,2),w(counter,3),'*')
hold off
az = az + 0.5;
view(az,el)
pause(0.01);
end
t(counter+1:end) = [];
w(counter+1,:) = [];
\ No newline at end of file
function [t,x] = RK4(t0,x0,h,F)
t = t0 + h;
k1 = F(t0,x0);
k2 = F(t0 + h/2, x0 + h.*k1./2);
k3 = F(t0 + h/2, x0 + h.*k2./2);
k4 = F(t0 + h, x0 + h.*k3);
x = x0 + (h/6).*(k1 + 2.*k2 + 2.*k3 + k4);
end
\ No newline at end of file
function varargout = Rotate(varargin)
% ROTATE MATLAB code for Rotate.fig
% ROTATE, by itself, creates a new ROTATE or raises the existing
% singleton*.
%
% H = ROTATE returns the handle to a new ROTATE or the handle to
% the existing singleton*.
%
% ROTATE('CALLBACK',hObject,eventData,handles,...) calls the local
% function named CALLBACK in ROTATE.M with the given input arguments.
%
% ROTATE('Property','Value',...) creates a new ROTATE or raises the
% existing singleton*. Starting from the left, property value pairs are
% applied to the GUI before Rotate_OpeningFcn gets called. An
% unrecognized property name or invalid value makes property application
% stop. All inputs are passed to Rotate_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 Rotate
% Last Modified by GUIDE v2.5 11-May-2016 14:43:26
% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name', mfilename, ...
'gui_Singleton', gui_Singleton, ...
'gui_OpeningFcn', @Rotate_OpeningFcn, ...
'gui_OutputFcn', @Rotate_OutputFcn, ...
'gui_LayoutFcn', [] , ...
'gui_Callback', []);
if nargin && ischar(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 Rotate is made visible.
function Rotate_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 Rotate (see VARARGIN)
% Choose default command line output for Rotate
handles.output = hObject;
handles.Az = 90;
handles.El = 0;
handles.Slider1Min = -180;
handles.Slider1Max = 180;
handles.Slider2Min = -90;
handles.Slider2Max = 90;
% Update handles structure
guidata(hObject, handles);
% UIWAIT makes Rotate wait for user response (see UIRESUME)
% uiwait(handles.figure1);
% --- Outputs from this function are returned to the command line.
function varargout = Rotate_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 slider movement.
function slider1_Callback(hObject, eventdata, handles)
% hObject handle to slider1 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
handles.Az = get(hObject,'Value'); %returns position of slider
get(hObject,'handles.Slider1Min');
get(hObject,'handles.Slider1Max'); %to determine range of slider
view(handles.Az,handles.El)
guidata(hObject,handles);
% --- Executes during object creation, after setting all properties.
function slider1_CreateFcn(hObject, eventdata, handles)
% hObject handle to slider1 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles empty - handles not created until after all CreateFcns called
% Hint: slider controls usually have a light gray background.
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
set(hObject,'BackgroundColor',[.9 .9 .9]);
end
% --- Executes on slider movement.
function slider2_Callback(hObject, eventdata, handles)
% hObject handle to slider2 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
handles.El = get(hObject,'Value'); %returns position of slider
get(hObject,'handles.Slider2Min');
get(hObject,'handles.Slider2Max'); %to determine range of slider
view(handles.Az,handles.El)
guidata(hObject,handles);
% --- Executes during object creation, after setting all properties.
function slider2_CreateFcn(hObject, eventdata, handles)
% hObject handle to slider2 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles empty - handles not created until after all CreateFcns called
% Hint: slider controls usually have a light gray background.
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
set(hObject,'BackgroundColor',[.9 .9 .9]);
end
% --- Executes during object creation, after setting all properties.
function axes1_CreateFcn(hObject, eventdata, handles)
% hObject handle to axes1 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles empty - handles not created until after all CreateFcns called
% Hint: place code in OpeningFcn to populate axes1
load('Sol.mat');
handles.w = w;
% Update handles structure
guidata(hObject, handles);
plot3(w(:,1),w(:,2),w(:,guide
3));
function varargout = Rotate2(varargin)
% ROTATE2 MATLAB code for Rotate2.fig
% ROTATE2, by itself, creates a new ROTATE2 or raises the existing
% singleton*.
%
% H = ROTATE2 returns the handle to a new ROTATE2 or the handle to
% the existing singleton*.
%
% ROTATE2('CALLBACK',hObject,eventData,handles,...) calls the local
% function named CALLBACK in ROTATE2.M with the given input arguments.
%
% ROTATE2('Property','Value',...) creates a new ROTATE2 or raises the
% existing singleton*. Starting from the left, property value pairs are
% applied to the GUI before Rotate2_OpeningFcn gets called. An
% unrecognized property name or invalid value makes property application
% stop. All inputs are passed to Rotate2_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 Rotate2
% Last Modified by GUIDE v2.5 12-May-2016 13:29:44
% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name', mfilename, ...
'gui_Singleton', gui_Singleton, ...
'gui_OpeningFcn', @Rotate2_OpeningFcn, ...
'gui_OutputFcn', @Rotate2_OutputFcn, ...
'gui_LayoutFcn', [] , ...
'gui_Callback', []);
if nargin && ischar(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 Rotate2 is made visible.
function Rotate2_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 Rotate2 (see VARARGIN)
% Choose default command line output for Rotate2
handles.output = hObject;
handles.Az = 90;
handles.El = 0;
handles.Slider1Min = -180;
handles.Slider1Max = 180;
handles.Slider2Min = -90;
handles.Slider2Max = 90;
handles.axes = [];
% Update handles structure
guidata(hObject, handles);
% UIWAIT makes Rotate2 wait for user response (see UIRESUME)
% uiwait(handles.figure1);
% --- Outputs from this function are returned to the command line.
function varargout = Rotate2_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 slider movement.
function slider1_Callback(hObject, eventdata, handles)
% hObject handle to slider1 (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,'Value') returns position of slider
% get(hObject,'Min') and get(hObject,'Max') to determine range of slider
handles.El = get(hObject,'Value');
get(hObject,'Min');
get(hObject,'Max'); %to determine range of slider
guidata(hObject, handles);
view(handles.Az,handles.El)
% --- Executes during object creation, after setting all properties.
function slider1_CreateFcn(hObject, eventdata, handles)
% hObject handle to slider1 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles empty - handles not created until after all CreateFcns called
% Hint: slider controls usually have a light gray background.
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
set(hObject,'BackgroundColor',[.9 .9 .9]);
end
% --- Executes during object creation, after setting all properties.
function axes1_CreateFcn(hObject, eventdata, handles)
% hObject handle to axes1 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles empty - handles not created until after all CreateFcns called
% Hint: place code in OpeningFcn to populate axes1
load('Sol.mat');
handles.w = w;
handles.axes = gca;
% Update handles structure
guidata(hObject, handles);
plot3(handles.axes,w(:,1),w(:,2),w(:,3));
view(0,90)
% --- Executes on slider movement.
function slider2_Callback(hObject, eventdata, handles)
% hObject handle to slider2 (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,'Value') returns position of slider
% get(hObject,'Min') and get(hObject,'Max') to determine range of slider
handles.Az = get(hObject,'Value');
get(hObject,'Min');
get(hObject,'Max'); %to determine range of slider
guidata(hObject, handles);
view(handles.Az,handles.El)
% --- Executes during object creation, after setting all properties.
function slider2_CreateFcn(hObject, eventdata, handles)
% hObject handle to slider2 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles empty - handles not created until after all CreateFcns called
% Hint: slider controls usually have a light gray background.
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
set(hObject,'BackgroundColor',[.9 .9 .9]);
end
{}
\ No newline at end of file
%!PS-Adobe-3.0 EPSF-3.0
%%Creator: MATLAB, The MathWorks, Inc. Version 8.1.0.604 (R2013a). Operating System: Linux 3.13.0-57-generic #95-Ubuntu SMP Fri Jun 19 09:28:15 UTC 2015 x86_64.
%%Title: /home/ggoedert/Dropbox/Mestrado GTG/2015.1/AnaliseNumedica/Lista5/Circle/CircleRK2C05.eps
%%CreationDate: 07/09/2015 17:35:25
%%DocumentNeededFonts: Helvetica-Bold
%%DocumentProcessColors: Cyan Magenta Yellow Black
%%LanguageLevel: 2
%%Pages: 1
%%BoundingBox: 81 227 529 564
%%EndComments
%%BeginProlog
% MathWorks dictionary
/MathWorks 160 dict begin
% definition operators
/bdef {bind def} bind def
/ldef {load def} bind def
/xdef {exch def} bdef
/xstore {exch store} bdef
% operator abbreviations
/c /clip ldef
/cc /concat ldef
/cp /closepath ldef
/gr /grestore ldef
/gs /gsave ldef
/mt /moveto ldef
/np /newpath ldef
/cm /currentmatrix ldef
/sm /setmatrix ldef
/rm /rmoveto ldef
/rl /rlineto ldef
/s {show newpath} bdef
/sc {setcmykcolor} bdef
/sr /setrgbcolor ldef
/sg /setgray ldef
/w /setlinewidth ldef
/j /setlinejoin ldef
/cap /setlinecap ldef
/rc {rectclip} bdef
/rf {rectfill} bdef
% page state control
/pgsv () def
/bpage {/pgsv save def} bdef
/epage {pgsv restore} bdef
/bplot /gsave ldef
/eplot {stroke grestore} bdef
% orientation switch
/portraitMode 0 def /landscapeMode 1 def /rotateMode 2 def
% coordinate system mappings
/dpi2point 0 def
% font control
/FontSize 0 def
/FMS {/FontSize xstore findfont [FontSize 0 0 FontSize neg 0 0]
makefont setfont} bdef
/reencode {exch dup where {pop load} {pop StandardEncoding} ifelse
exch dup 3 1 roll findfont dup length dict begin
{ 1 index /FID ne {def}{pop pop} ifelse } forall
/Encoding exch def currentdict end definefont pop} bdef
/isroman {findfont /CharStrings get /Agrave known} bdef
/FMSR {3 1 roll 1 index dup isroman {reencode} {pop pop} ifelse
exch FMS} bdef
/csm {1 dpi2point div -1 dpi2point div scale neg translate
dup landscapeMode eq {pop -90 rotate}
{rotateMode eq {90 rotate} if} ifelse} bdef
% line types: solid, dotted, dashed, dotdash
/SO { [] 0 setdash } bdef
/DO { [.5 dpi2point mul 4 dpi2point mul] 0 setdash } bdef
/DA { [6 dpi2point mul] 0 setdash } bdef
/DD { [.5 dpi2point mul 4 dpi2point mul 6 dpi2point mul 4
dpi2point mul] 0 setdash } bdef
% macros for lines and objects
/L {lineto stroke} bdef
/MP {3 1 roll moveto 1 sub {rlineto} repeat} bdef
/AP {{rlineto} repeat} bdef
/PDlw -1 def
/W {/PDlw currentlinewidth def setlinewidth} def
/PP {closepath eofill} bdef
/DP {closepath stroke} bdef
/MR {4 -2 roll moveto dup 0 exch rlineto exch 0 rlineto
neg 0 exch rlineto closepath} bdef
/FR {MR stroke} bdef
/PR {MR fill} bdef
/L1i {{currentfile picstr readhexstring pop} image} bdef
/tMatrix matrix def
/MakeOval {newpath tMatrix currentmatrix pop translate scale
0 0 1 0 360 arc tMatrix setmatrix} bdef
/FO {MakeOval stroke} bdef
/PO {MakeOval fill} bdef
/PD {currentlinewidth 2 div 0 360 arc fill
PDlw -1 eq not {PDlw w /PDlw -1 def} if} def
/FA {newpath tMatrix currentmatrix pop translate scale
0 0 1 5 -2 roll arc tMatrix setmatrix stroke} bdef
/PA {newpath tMatrix currentmatrix pop translate 0 0 moveto scale
0 0 1 5 -2 roll arc closepath tMatrix setmatrix fill} bdef
/FAn {newpath tMatrix currentmatrix pop translate scale
0 0 1 5 -2 roll arcn tMatrix setmatrix stroke} bdef
/PAn {newpath tMatrix currentmatrix pop translate 0 0 moveto scale
0 0 1 5 -2 roll arcn closepath tMatrix setmatrix fill} bdef
/vradius 0 def /hradius 0 def /lry 0 def
/lrx 0 def /uly 0 def /ulx 0 def /rad 0 def
/MRR {/vradius xdef /hradius xdef /lry xdef /lrx xdef /uly xdef
/ulx xdef newpath tMatrix currentmatrix pop ulx hradius add uly
vradius add translate hradius vradius scale 0 0 1 180 270 arc
tMatrix setmatrix lrx hradius sub uly vradius add translate
hradius vradius scale 0 0 1 270 360 arc tMatrix setmatrix
lrx hradius sub lry vradius sub translate hradius vradius scale
0 0 1 0 90 arc tMatrix setmatrix ulx hradius add lry vradius sub
translate hradius vradius scale 0 0 1 90 180 arc tMatrix setmatrix
closepath} bdef
/FRR {MRR stroke } bdef
/PRR {MRR fill } bdef
/MlrRR {/lry xdef /lrx xdef /uly xdef /ulx xdef /rad lry uly sub 2 div def
newpath tMatrix currentmatrix pop ulx rad add uly rad add translate
rad rad scale 0 0 1 90 270 arc tMatrix setmatrix lrx rad sub lry rad
sub translate rad rad scale 0 0 1 270 90 arc tMatrix setmatrix
closepath} bdef
/FlrRR {MlrRR stroke } bdef
/PlrRR {MlrRR fill } bdef
/MtbRR {/lry xdef /lrx xdef /uly xdef /ulx xdef /rad lrx ulx sub 2 div def
newpath tMatrix currentmatrix pop ulx rad add uly rad add translate
rad rad scale 0 0 1 180 360 arc tMatrix setmatrix lrx rad sub lry rad
sub translate rad rad scale 0 0 1 0 180 arc tMatrix setmatrix
closepath} bdef
/FtbRR {MtbRR stroke } bdef
/PtbRR {MtbRR fill } bdef
/stri 6 array def /dtri 6 array def
/smat 6 array def /dmat 6 array def
/tmat1 6 array def /tmat2 6 array def /dif 3 array def
/asub {/ind2 exch def /ind1 exch def dup dup
ind1 get exch ind2 get sub exch } bdef
/tri_to_matrix {
2 0 asub 3 1 asub 4 0 asub 5 1 asub
dup 0 get exch 1 get 7 -1 roll astore } bdef
/compute_transform {
dmat dtri tri_to_matrix tmat1 invertmatrix
smat stri tri_to_matrix tmat2 concatmatrix } bdef
/ds {stri astore pop} bdef
/dt {dtri astore pop} bdef
/db {2 copy /cols xdef /rows xdef mul dup 3 mul string
currentfile
3 index 0 eq {/ASCIIHexDecode filter}