(original) (raw)

gnimatlab/ 40755 12122 12120 0 7553235722 11463 5ustar hairermathgnimatlab/coeff_comp.m100755 12122 12120 10006 7551506025 14051 0ustar hairermathfunction gamma = coeff_comp(name) %COEFF_COMP Set coefficients for the GNI_COMP integrator % Syntax: GAMMA = COEFF_COMP('NAME'); % % See also GNI_COMP switch name case '43', gamma = ones(3,1); qr2 = 2^(1/3); gamma(1) = 1/(2-qr2); gamma(2) = -qr2*gamma(1); gamma(3) = gamma(1); case '45', gamma = ones(5,1); qr4 = 4^(1/3); gamma(1) = 1/(4-qr4); gamma(2) = gamma(1); gamma(3) = -qr4*gamma(1); gamma(4) = gamma(1); gamma(5) = gamma(1); case '67', gamma = ones(7,1); gamma(1)= 0.78451361047755726381949763; gamma(2)= 0.23557321335935813368479318; gamma(3)=-1.17767998417887100694641568; gamma(4)= 1.31518632068391121888424973; gamma(5)=-1.17767998417887100694641568; gamma(6)= 0.23557321335935813368479318; gamma(7)= 0.78451361047755726381949763; case '69', gamma = ones(9,1); gamma(1)= 0.39215846226558135393725811; gamma(2)= 0.33260288849639382848984234; gamma(3)=-0.70624607359230088534941023; gamma(4)= 0.08221272729243641010472095; gamma(5)= 0.79854399107577858563517766; gamma(6)= 0.08221272729243641010472095; gamma(7)=-0.70624607359230088534941023; gamma(8)= 0.33260288849639382848984234; gamma(9)= 0.39215846226558135393725811; case '815', gamma = ones(8,1); gamma(1)= 0.74167036435061295344822780; gamma(2)=-0.40910082580003159399730010; gamma(3)= 0.19075471029623837995387626; gamma(4)=-0.57386247111608226665638773; gamma(5)= 0.29906418130365592384446354; gamma(6)= 0.33462491824529818378495798; gamma(7)= 0.31529309239676659663205666; gamma(8)=-0.79688793935291635401978884; gamma(9)= 0.31529309239676659663205666; gamma(10)= 0.33462491824529818378495798; gamma(11)= 0.29906418130365592384446354; gamma(12)=-0.57386247111608226665638773; gamma(13)= 0.19075471029623837995387626; gamma(14)=-0.40910082580003159399730010; gamma(15)= 0.74167036435061295344822780; case '817', gamma = ones(17,1); gamma( 1)= 0.13020248308889008087881763; gamma( 2)= 0.56116298177510838456196441; gamma( 3)=-0.38947496264484728640807860; gamma( 4)= 0.15884190655515560089621075; gamma( 5)=-0.39590389413323757733623154; gamma( 6)= 0.18453964097831570709183254; gamma( 7)= 0.25837438768632204729397911; gamma( 8)= 0.29501172360931029887096624; gamma( 9)=-0.60550853383003451169892108; gamma(10)= 0.29501172360931029887096624; gamma(11)= 0.25837438768632204729397911; gamma(12)= 0.18453964097831570709183254; gamma(13)=-0.39590389413323757733623154; gamma(14)= 0.15884190655515560089621075; gamma(15)=-0.38947496264484728640807860; gamma(16)= 0.56116298177510838456196441; gamma(17)= 0.13020248308889008087881763; case '1033', gamma = ones(33,1); gamma( 1)= 0.09040619368607278492161150; gamma( 2)= 0.53591815953030120213784983; gamma( 3)= 0.35123257547493978187517736; gamma( 4)=-0.31116802097815835426086544; gamma( 5)=-0.52556314194263510431065549; gamma( 6)= 0.14447909410225247647345695; gamma( 7)= 0.02983588609748235818064083; gamma( 8)= 0.17786179923739805133592238; gamma( 9)= 0.09826906939341637652532377; gamma(10)= 0.46179986210411860873242126; gamma(11)=-0.33377845599881851314531820; gamma(12)= 0.07095684836524793621031152; gamma(13)= 0.23666960070126868771909819; gamma(14)=-0.49725977950660985445028388; gamma(15)=-0.30399616617237257346546356; gamma(16)= 0.05246957188100069574521612; gamma(17)= 0.44373380805019087955111365; gamma(18)= 0.05246957188100069574521612; gamma(19)=-0.30399616617237257346546356; gamma(20)=-0.49725977950660985445028388; gamma(21)= 0.23666960070126868771909819; gamma(22)= 0.07095684836524793621031152; gamma(23)=-0.33377845599881851314531820; gamma(24)= 0.46179986210411860873242126; gamma(25)= 0.09826906939341637652532377; gamma(26)= 0.17786179923739805133592238; gamma(27)= 0.02983588609748235818064083; gamma(28)= 0.14447909410225247647345695; gamma(29)=-0.52556314194263510431065549; gamma(30)=-0.31116802097815835426086544; gamma(31)= 0.35123257547493978187517736; gamma(32)= 0.53591815953030120213784983; gamma(33)= 0.09040619368607278492161150; case '21', gamma = 1; otherwise, gamma = 1; error(['The method ''' lower(name) ''' is not implemented.']); end gnimatlab/coeff_irk2.m100755 12122 12120 20520 7544612461 13770 0ustar hairermathfunction [ns,nm,C,B,BC,AA,E,SM,AM] = coeff_irk2(method) nm = 3; switch method case 'G4', ns = 2; C = zeros(ns,1); B = zeros(ns,1); BC = zeros(ns,1); AA = zeros(ns,ns); E = zeros(ns,ns+nm); SM = zeros(nm,1); AM = zeros(ns + nm,1); C(1) = 0.21132486540518711775E+00; C(2) = 0.78867513459481288225E+00; B(1) = 0.50000000000000000000E+00; B(2) = 0.50000000000000000000E+00; BC(1) = 0.39433756729740644113E+00; BC(2) = 0.10566243270259355887E+00; AA(1,1)= 0.41666666666666666667E-01; AA(1,2)=-0.19337567297406441127E-01; AA(2,1)= 0.26933756729740644113E+00; AA(2,2)= 0.41666666666666666667E-01; E(1,1)=-0.28457905077110526160E-02; E(1,2)=-0.63850024471784160410E-01; E(1,3)= 0.48526095198694517563E-02; E(1,4)= 0.11305688530429939012E+00; E(1,5)=-0.28884580475413403312E-01; E(2,1)= 0.41122751744511433137E-01; E(2,2)=-0.18654814888622834132E+00; E(2,3)=-0.18110185277445209332E-01; E(2,4)= 0.36674109449368040786E+00; E(2,5)= 0.10779872188955481745E+00; SM(1)= 0.00000000000000000000E+00; SM(2)= 0.10000000000000000000E+01; SM(3)= 0.16000000000000000000E+01; AM(1)= 0.25279583039343438291E+02; AM(2)=-0.86907830393434382912E+01; AM(3)=-0.80640000000000000000E+00; AM(4)= 0.29184000000000000000E+01; AM(5)= 0.00000000000000000000E+00; case 'G8', ns = 4; C = zeros(ns,1); B = zeros(ns,1); BC = zeros(ns,1); AA = zeros(ns,ns); E = zeros(ns,ns+nm); SM = zeros(nm,1); AM = zeros(ns + nm,1); C(1)= 0.69431844202973712388E-01; C(2)= 0.33000947820757186760E+00; C(3)= 0.66999052179242813240E+00; C(4)= 0.93056815579702628761E+00; B(1)= 0.17392742256872692869E+00; B(2)= 0.32607257743127307131E+00; B(3)= 0.32607257743127307131E+00; B(4)= 0.17392742256872692869E+00; BC(1)= 0.16185132086231030665E+00; BC(2)= 0.21846553629538057030E+00; BC(3)= 0.10760704113589250101E+00; BC(4)= 0.12076101706416622036E-01; AA(1,1)= 0.40381914508467311298E-02; AA(1,2)=-0.32958609449446961650E-02; AA(1,3)= 0.26447829520668538006E-02; AA(1,4)=-0.97672296325588161023E-03; AA(2,1)= 0.43563580902396261254E-01; AA(2,2)= 0.13818951406296126013E-01; AA(2,3)=-0.43401341944349953440E-02; AA(2,4)= 0.14107297391595337720E-02; AA(3,1)= 0.10586435263357640763E+00; AA(3,2)= 0.10651836096505307395E+00; AA(3,3)= 0.13818951406296126013E-01; AA(3,4)=-0.17580153590805494993E-02; AA(4,1)= 0.14879849619263780300E+00; AA(4,2)= 0.19847049885237718995E+00; AA(4,3)= 0.81671359795877570687E-01; AA(4,4)= 0.40381914508467311298E-02; E(1,1)=-0.21272768296134340207E-01; E(1,2)= 0.11059138674756969912E-01; E(1,3)= 0.38999255049973564023E-02; E(1,4)=-0.43986226789008967612E-01; E(1,5)= 0.13581590305438849621E-01; E(1,6)= 0.39922421675314269059E-01; E(1,7)=-0.79369058065113002021E-03; E(2,1)=-0.75671119283734809953E-02; E(2,2)= 0.10209394000843457002E-01; E(2,3)=-0.12880197817980892596E-01; E(2,4)=-0.56381316813776501277E-01; E(2,5)= 0.37440782682669799960E-02; E(2,6)= 0.11522469441011273193E+00; E(2,7)= 0.21035877343246316334E-02; E(3,1)=-0.39890571772473709759E+00; E(3,2)= 0.26819725655216894347E+00; E(3,3)=-0.82551711648854471247E-01; E(3,4)=-0.85516559106259630212E+00; E(3,5)= 0.24433810515772642570E+00; E(3,6)= 0.10234155624049009806E+01; E(3,7)= 0.25115745967236579242E-01; E(4,1)=-0.40964796048052939224E+00; E(4,2)= 0.29949323098224574487E+00; E(4,3)=-0.13867460566101912494E+00; E(4,4)=-0.98859300714628940382E+00; E(4,5)= 0.24671351779481625627E+00; E(4,6)= 0.12912760231350872304E+01; E(4,7)= 0.13241134766742798418E+00; SM(1)= 0.00000000000000000000E+00; SM(2)= 0.10000000000000000000E+01; SM(3)= 0.16500000000000000000E+01; AM(1)= 0.10806374869244001787E+04; AM(2)=-0.66008818661284690206E+03; AM(3)= 0.61810154357557529566E+03; AM(4)=-0.31341427826212857229E+03; AM(5)=-0.10187174765625000000E+02; AM(6)= 0.31173050390625000000E+02; AM(7)= 0.00000000000000000000E+00; case 'G12', ns = 6; C = zeros(ns,1); B = zeros(ns,1); BC = zeros(ns,1); AA = zeros(ns,ns); E = zeros(ns,ns+nm); SM = zeros(nm,1); AM = zeros(ns + nm,1); C(1)= 0.33765242898423986094E-01; C(2)= 0.16939530676686774317E+00; C(3)= 0.38069040695840154568E+00; C(4)= 0.61930959304159845432E+00; C(5)= 0.83060469323313225683E+00; C(6)= 0.96623475710157601391E+00; B(1)= 0.85662246189585172520E-01; B(2)= 0.18038078652406930378E+00; B(3)= 0.23395696728634552369E+00; B(4)= 0.23395696728634552369E+00; B(5)= 0.18038078652406930378E+00; B(6)= 0.85662246189585172520E-01; BC(1)= 0.82769839639769234611E-01; BC(2)= 0.14982512785597570103E+00; BC(3)= 0.14489179419935320895E+00; BC(4)= 0.89065173086992314743E-01; BC(5)= 0.30555658668093602753E-01; BC(6)= 0.28924065498159379092E-02; AA(1,1)= 0.90625420195651151857E-03; AA(1,2)=-0.72859711612531400024E-03; AA(1,3)= 0.79102695861167691135E-03; AA(1,4)=-0.70675390218535384182E-03; AA(1,5)= 0.45647714224056921122E-03; AA(1,6)=-0.14836147050330408643E-03; AA(2,1)= 0.11272367531794365387E-01; AA(2,2)= 0.39083482447840698486E-02; AA(2,3)=-0.14724868010943911900E-02; AA(2,4)= 0.10992669056588431310E-02; AA(2,5)=-0.67689040729401428165E-03; AA(2,6)= 0.21677950347174141516E-03; AA(3,1)= 0.30008019623627547434E-01; AA(3,2)= 0.36978289259468146662E-01; AA(3,3)= 0.65490339168957822692E-02; AA(3,4)=-0.16615098173008262274E-02; AA(3,5)= 0.84753461862041607649E-03; AA(3,6)=-0.25877462623437421721E-03; AA(4,1)= 0.49900269650650898941E-01; AA(4,2)= 0.82003427445271620462E-01; AA(4,3)= 0.54165111295060067982E-01; AA(4,4)= 0.65490339168957822692E-02; AA(4,5)=-0.11352871017627472322E-02; AA(4,6)= 0.28963081055952389031E-03; AA(5,1)= 0.68475836671617248304E-01; AA(5,2)= 0.11859257878058808400E+00; AA(5,3)= 0.10635984886129551097E+00; AA(5,4)= 0.47961474042181382443E-01; AA(5,5)= 0.39083482447840698486E-02; AA(5,6)=-0.34600839001342442657E-03; AA(6,1)= 0.79729071619449992615E-01; AA(6,2)= 0.14419100392702230613E+00; AA(6,3)= 0.13628542646896576408E+00; AA(6,4)= 0.81956586217401900627E-01; AA(6,5)= 0.23736460480774324642E-01; AA(6,6)= 0.90625420195651151857E-03; E(1,1)=-0.16761132335280609813E-01; E(1,2)= 0.10201050166615899799E-01; E(1,3)=-0.58593121685075943100E-02; E(1,4)=-0.11907383391366998251E-03; E(1,5)= 0.10615611118132982241E-01; E(1,6)=-0.30692054230989138447E-01; E(1,7)= 0.10615182045216224925E-01; E(1,8)= 0.22586707045496892369E-01; E(1,9)=-0.16931992776201068110E-04; E(2,1)= 0.10671755276327262128E-01; E(2,2)=-0.51098203653251450913E-02; E(2,3)= 0.16062647299186369205E-03; E(2,4)= 0.64818802653621866868E-02; E(2,5)=-0.12132386914873895089E-01; E(2,6)=-0.99709737725909584834E-02; E(2,7)=-0.70287093442894942752E-02; E(2,8)= 0.31243249755879001843E-01; E(2,9)= 0.31763603839792897936E-04; E(3,1)=-0.40875203230945019464E+00; E(3,2)= 0.28214948905763253599E+00; E(3,3)=-0.22612660499718519054E+00; E(3,4)= 0.13640993962985420478E+00; E(3,5)= 0.15888529591697266925E+00; E(3,6)=-0.11667863471317749710E+01; E(3,7)= 0.25224964119340060668E+00; E(3,8)= 0.10440940643938620983E+01; E(3,9)= 0.33914722176493324285E-03; E(4,1)=-0.29437531285359759661E+01; E(4,2)= 0.20017220470127690267E+01; E(4,3)=-0.15383035791443948798E+01; E(4,4)= 0.78114323215109899716E+00; E(4,5)= 0.13930345104184182146E+01; E(4,6)=-0.75958281612589849630E+01; E(4,7)= 0.18220129530415584951E+01; E(4,8)= 0.62663163493155487560E+01; E(4,9)= 0.54279630166374655267E-02; E(5,1)=-0.79572842006457093076E+01; E(5,2)= 0.53527892762707449170E+01; E(5,3)=-0.40049139768467199697E+01; E(5,4)= 0.18326058141135591515E+01; E(5,5)= 0.39753886181058367500E+01; E(5,6)=-0.19423696478604790213E+02; E(5,7)= 0.49362128400107292627E+01; E(5,8)= 0.15601708062381928560E+02; E(5,9)= 0.32142123424873719685E-01; E(6,1)=-0.78463118056075171475E+01; E(6,2)= 0.53580869574441241664E+01; E(6,3)=-0.41476905275607763365E+01; E(6,4)= 0.21275912797813913113E+01; E(6,5)= 0.37642416878253538582E+01; E(6,6)=-0.20329681631523484613E+02; E(6,7)= 0.48515418060343387549E+01; E(6,8)= 0.16604467346259915039E+02; E(6,9)= 0.84559690262225766975E-01; SM(1)= 0.00000000000000000000E+00; SM(2)= 0.10000000000000000000E+01; SM(3)= 0.17500000000000000000E+01; AM(1)= 0.58080578375796358720E+05; AM(2)=-0.33214989339522861968E+05; AM(3)= 0.28376088288312020853E+05; AM(4)=-0.27923430684614999462E+05; AM(5)= 0.29743005589491042677E+05; AM(6)=-0.15525927919158826444E+05; AM(7)=-0.27700591278076171875E+03; AM(8)= 0.73086943817138671875E+03; AM(9)= 0.00000000000000000000E+00; otherwise error(['The required method "' method '" is not implemented']); end gnimatlab/coeff_lmm2.m100755 12122 12120 557 7544612431 13735 0ustar hairermathfunction [A,B,C] = coeff_lmm2(method) switch method case '801', A = [1;0;1;0.5]; B = [17671;-23622;61449;-25258]/12096; case '802', A = [1;2;3;1.75]; B = [192481;6582;816783;-78406]/120960; case '803', A = [1;1;1;0.5]; B = [13207;-8934;42873;-16906]/8640; otherwise error(['The required method "' method '" is not implemented']); end C = [672;-168;32;-3]/840; gnimatlab/Contents.m100644 12122 12120 1343 7553235514 13533 0ustar hairermath% GniMatlab -- Geometric Numerical Integration codes for Matlab. % % GNI Integrators. % gni_irk2 - Implicit Runge-Kutta methods. % gni_lmm2 - Symmetric multistep methods. % gni_comp - Generic composition methods. % % GNI Option handling. % gniset - Create/alter GNI OPTIONS structure. % gniget - Get GNI OPTIONS parameters. % % Examples of GNI problem files. % henon - Henon-Heiles problem. % kepler - Kepler problem. % twobodysphere - Two-body problem on a sphere. % % Examples of custom output functions. % phaseplot - Two-dimensional plot in phase space. % sphereplot - Three-dimensional plot on a sphere. % Copyright (c) 2002 by Ernst Hairer & Martin Hairer. gnimatlab/gni_comp.m100644 12122 12120 15756 7553206013 13556 0ustar hairermathfunction [tout,qout,pout,teout,qeout,peout,ieout] = gni_comp(odefile,tspan,y0,options,varargin) %GNI_COMP Solves ordinary differential equations with a composition method. % The calling sequence is identical to GNI_IRK2, see HELP GNI_IRK2 for details. % % The standard call uses a Stoermer/Verlet scheme and applies to second order % differential equations, but the basic mathod can be changed with the % 'PartialFlow' option, so that general problems can be solved. % % Details are to be found in the article "GniCodes - Matlab Programs for Geometric % Numerical Integration", available from http://www.unige.ch/math/folks/hairer/ % under item software % % Copyright by Ernst Hairer and Martin Hairer, 10.10.2002. true = logical(1); false = ~true; nsteps = 0; % stats nfailed = 0; % stats nfevals = 0; % stats npds = 0; % stats ndecomps = 0; % stats nsolves = 0; % stats if nargin == 0 error('Not enough input arguments. See GNI_COMP.'); elseif ~isstr(odefile) & ~isa(odefile, 'inline') error('First argument must be a single-quoted string. See GNI_COMP.'); end if nargin == 1 tspan = []; y0 = []; options = []; elseif nargin == 2 y0 = []; options = []; elseif nargin == 3 options = []; elseif ~isempty(options) & ~isa(options,'struct') if (length(tspan) == 1) & (length(y0) == 1) & (min(size(options)) == 1) msg = sprintf('Use gni_comp(''%s'',tspan,y0,...) instead.',odefile); error(['Obsolete syntax. ' msg]); else error('Correct syntax is gni_comp(''odefile'',tspan,y0,options).'); end end if isstr(odefile) & (exist(odefile) == 4) % a SIMULINK model error('Use the SIM command directly for SIMULINK models.'); return; end % Get default tspan and y0 from odefile if none are specified. if isempty(tspan) | isempty(y0) if (nargout(odefile) < 3) & (nargout(odefile) ~= -1) msg = sprintf('Use gni_comp(''%s'',tspan,y0,...) instead.',odefile); error(['No default parameters in ' upper(odefile) '. ' msg]); end [def_tspan,def_y0,def_options] = feval(odefile,[],[],'init',varargin{:}); if isempty(tspan) tspan = def_tspan; end if isempty(y0) y0 = def_y0; end if isempty(options) options = def_options; else options = gniset(def_options,options); end end % Test that tspan is internally consistent. if size(y0,2) ~= 1 y0 = y0'; end ntspan = length(tspan); if ntspan > 2 error('The timespan must consist of one or two numbers.'); end if ntspan == 1 tspan = [0;tspan]; t0 = 0; else t0 = tspan(1); end tfinal = tspan(2); if t0 >= tfinal error('The final time must greater than the starting time.'); end outstep = gniget(options,'OutputSteps',1); if (size(outstep) ~= [1 1]) error('The option ''OutputSteps'' must contain a single number'); end outflag = 0; if (outstep > 0) outflag = 1; else outstep = 1; end t = t0; if (mod(length(y0),2) ~= 0) error('The initial value must have 2*N components, where N is the dimension of the system.'); end neq = length(y0)/2; Q = y0(1:neq); P = y0((neq+1):(2*neq)); locateevents = strcmp(gniget(options,'Events','off'),'on'); teout = []; qeout = []; peout = []; ieout = []; if locateevents [V,terminal,direction] = feval(odefile,t,[Q;P],'events',varargin{:}); end if nargout > 0 outfun = gniget(options,'OutputFcn'); else outfun = gniget(options,'OutputFcn','odeplot'); end if isempty(outfun) haveoutfun = false; else haveoutfun = true; outputs = gniget(options,'OutputSel',1:(2*neq)); end h = gniget(options,'StepSize'); if (isempty(h)) nsteps = gniget(options,'NumSteps'); if (isempty(nsteps)) h = 1e-2; fprintf('Warning: No initial step size provided, using h=1e-2 instead.\n'); else h = (tfinal - t0)/nsteps; end end if ((~isa(h,'double')) | (size(h) ~= [1 1])) error('The option ''StepSize'' must contain a single number'); end % Allocate memory if we're generating output. chunk = 1; if nargout > 0 chunk = max(ceil(128 / neq),1); tout = zeros(chunk,1); pout = zeros(chunk,neq); qout = zeros(chunk,neq); nout = 1; tout(nout) = t; qout(nout,:) = Q.'; pout(nout,:) = P.'; end % Initialize the output function. if haveoutfun QP = [Q;P]; feval(outfun,[t tfinal],QP(outputs),'init'); end % Initialize the IRK method basic = gniget(options,'PartialFlow','stverl'); if isempty(varargin) args = {}; else args = [{[]} varargin]; end feval(basic,t,0,0,odefile,0,0,0,false,false,'init',args{:}); method = gniget(options,'Method','817'); gamma = coeff_comp(method); ng = size(gamma); ng = max(ng); % THE MAIN LOOP done = false; outpoint = false; addpoint = false; HG = h*gamma; HA = HG(1) / 2; HGP = (HG(1:(ng-1)) + HG(2:ng))/2; HC = HG(ng) / 2; HGP = [HA;HGP;HC]; HFL = HA + HC; steppos = 0; first = true; last = false; nsteps = round((tfinal - t0)/h); for iter = 1:nsteps tnew = t+h; addpoint = false; outpoint = outpoint + 1; if (outpoint >= outstep) outpoint = 0; addpoint = true; end if locateevents QOld = Q; POld = P; VOld = V; end % Start of actual step ynew = y + h*... for i = 1:(ng-1) [P,Q] = feval(basic,t,P,Q,odefile,HGP(i),HG(i),HGP(i+1),first,last,[],args{:}); first = false; end if (haveoutfun | addpoint | locateevents) last = true; [P,Q] = feval(basic,t,P,Q,odefile,HGP(ng),HG(ng),HGP(ng+1),first,last,[],args{:}); first = true; last = false; else if done last = true; [P,Q] = feval(basic,t,P,Q,odefile,HGP(ng),HG(ng),HGP(ng+1),first,last,[],args{:}); else [P,Q] = feval(basic,t,P,Q,odefile,HGP(ng),HG(ng),HFL,first,last,[],args{:}); first = false; end end steppos = steppos+1; t = t0+steppos*h; if nargout > 0 if addpoint & outflag == 1 nout = nout + 1; if nout > length(tout) tout = [tout; zeros(chunk,1)]; pout = [pout; zeros(chunk,neq)]; qout = [qout; zeros(chunk,neq)]; end tout(nout) = t; qout(nout,:) = Q.'; pout(nout,:) = P.'; end if haveoutfun & addpoint QP = [Q;P]; if feval(outfun,t,QP(outputs),'') == 1 return; end end elseif haveoutfun & addpoint QP = [Q;P]; if feval(outfun,t,QP(outputs),'') == 1 return; end end if locateevents V = feval(odefile,t,[Q;P],'events',varargin{:}); [nteo,nieo,nqeo,npeo,stop] = gnievents(odefile,VOld,V,POld,P,QOld,Q,[t-h,t],direction,terminal,varargin{:}); teout = [teout;nteo]; ieout = [ieout;nieo]; qeout = [qeout;nqeo]; peout = [peout;npeo]; if stop term = true; break; end end % If there were no failures compute a new h. % Advance the integration one step. end feval(basic,t,0,0,odefile,0,0,0,false,false,'done',args{:}); if (nargout > 0 & outflag == 0) nout = nout + 1; if nout > length(tout) tout = [tout; zeros(1,1)]; pout = [pout; zeros(1,neq)]; qout = [qout; zeros(1,neq)]; end tout(nout) = t; qout(nout,:) = Q.'; pout(nout,:) = P.'; end if haveoutfun feval(outfun,[],[],'done'); end if nargout > 0 tout = tout(1:nout); qout = qout(1:nout,:); pout = pout(1:nout,:); end gnimatlab/gni_irk2.m100644 12122 12120 22357 7553205677 13501 0ustar hairermathfunction [tout,qout,pout,teout,qeout,peout,ieout,term] = gni_irk2(odefile,tspan,y0,options,varargin) %GNI_IRK2 Solves second-order differential equations with an implicit Runge-Kutta method. % [T,Q,P] = GNI_IRK2('F',TSPAN,[Q0;P0]) with TSPAN = [T0 TFINAL] integrates the % system of differential equations q'' = F(t,q) from time T0 to TFINAL with % initial conditions q = Q0, q' = P0. 'F' is a string containing the name of % a GNI problem file. Function F(t,q) must return a column vector. Each row in % solution arrays Q, P corresponds to a time returned in column vector T. % % [T,Q,P] = GNI_IRK2('F',TSPAN,[Q0;P0],OPTIONS) solves as above with default % integration parameters replaced by values in OPTIONS, an argument % created with the GNISET function. See GNISET for details. Additional % parameters can be passed to the GNI problem file as in the ODE suite. % % It is also possible to specify TSPAN, Y0 and/or OPTIONS in the GNI problem % file. In this case, the corresponding arguments should be empty. % % As an example, the command % % gni_irk2('kepler',[],[],[],0.5); % % solves the two-body Kepler problem with default parameters defined in % kepler.m and eccentricity equal to 0.5. % % [T,Q,P,TE,PE,QE,IE] = GNI_IRK2('F',TSPAN,[Q0;P0],OPTIONS) with the Events % property in OPTIONS set to 'on', solves as above while also locating zero % crossings of event function(s) defined in the GNI problem file. As an example, % the sequence of commands % % [T,Q,P,TE,QE,PE] = gni_irk2('henon'); % plot3(QE(:,2),PE(:,1),PE(:,2),'o'); % % draws the Poincare section of the standard Henon-Heiles problem. % % GNI_IRK2 implements symmetric symplectic Gauss methods of order 4, 8, and 12. % % Details are to be found in the article "GniCodes - Matlab Programs for Geometric % Numerical Integration", available from http://www.unige.ch/math/folks/hairer/ % under item software % % Copyright by Ernst Hairer and Martin Hairer, 10.10.2002. true = logical(1); false = ~true; term = false; nstps = 0; if nargin == 0 error('Not enough input arguments. See GNI_IRK2.'); elseif ~isstr(odefile) & ~isa(odefile, 'inline') error('First argument must be a single-quoted string. See GNI_IRK2.'); end if nargin == 1 tspan = []; y0 = []; options = []; elseif nargin == 2 y0 = []; options = []; elseif nargin == 3 options = []; elseif ~isempty(options) & ~isa(options,'struct') if (length(tspan) == 1) & (length(y0) == 1) & (min(size(options)) == 1) msg = sprintf('Use gni_irk2(''%s'',tspan,y0,...) instead.',odefile); error(['Obsolete syntax. ' msg]); else error('Correct syntax is gni_irk2(''odefile'',tspan,y0,options).'); end end if isstr(odefile) & (exist(odefile) == 4) % a SIMULINK model error('Use the SIM command directly for SIMULINK models.'); return; end % Get default tspan and y0 from odefile if none are specified. if isempty(tspan) | isempty(y0) if (nargout(odefile) < 3) & (nargout(odefile) ~= -1) msg = sprintf('Use gni_irk2(''%s'',tspan,y0,...) instead.',odefile); error(['No default parameters in ' upper(odefile) '. ' msg]); end [def_tspan,def_y0,def_options] = feval(odefile,[],[],'init',varargin{:}); if isempty(tspan) tspan = def_tspan; end if isempty(y0) y0 = def_y0; end if isempty(options) options = def_options; else options = gniset(def_options,options); end end % Test that tspan is internally consistent. if size(y0,2) ~= 1 y0 = y0'; end ntspan = length(tspan); if ntspan > 2 error('The timespan must consist of one or two numbers.'); end if ntspan == 1 tspan = [0;tspan]; t0 = 0; else t0 = tspan(1); end tfinal = tspan(2); if t0 >= tfinal error('The final time must greater than the starting time.'); end outstep = gniget(options,'OutputSteps',1); if (size(outstep) ~= [1 1]) error('The option ''OutputSteps'' must contain a single number'); end outflag = 0; if (outstep > 0) outflag = 1; else outstep = 1; end t = t0; if (mod(length(y0),2) ~= 0) error('The initial value must have 2*N components, where N is the dimension of the system.'); end neq = length(y0)/2; Q = y0(1:neq); P = y0((neq+1):(2*neq)); if nargout > 0 outfun = gniget(options,'OutputFcn'); else outfun = gniget(options,'OutputFcn','odeplot'); end if isempty(outfun) haveoutfun = false; else haveoutfun = true; outputs = gniget(options,'OutputSel',1:(2*neq)); end canvector = strcmp(gniget(options,'Vectorized','off'),'on'); locateevents = strcmp(gniget(options,'Events','off'),'on'); teout = []; qeout = []; peout = []; ieout = []; if locateevents [V,terminal,direction] = feval(odefile,t,[Q;P],'events',varargin{:}); end maxiter = gniget(options,'MaxIter',50); h = gniget(options,'StepSize'); if (isempty(h)) nsteps = gniget(options,'NumSteps'); if (isempty(nsteps)) h = 1e-2; fprintf('Warning: No initial step size provided, using h=1e-2 instead.\n'); else h = (tfinal - t0)/nsteps; end end if ((~isa(h,'double')) | (size(h) ~= [1 1])) error('The option ''StepSize'' must contain a single number'); end % Allocate memory if we're generating output. chunk = 1; if nargout > 0 chunk = max(ceil(128 / neq),1); tout = zeros(chunk,1); pout = zeros(chunk,neq); qout = zeros(chunk,neq); nout = 1; tout(nout) = t; qout(nout,:) = Q.'; pout(nout,:) = P.'; end if isempty(varargin) args = {}; else args = [{[]} varargin]; end FS = feval(odefile,t,Q,args{:}); [m,n] = size(FS); if n > 1 error([upper(odefile) ' must return a column vector.']) elseif m ~= neq msg = sprintf('an initial condition vector of length 2*%d.',m); error(['Solving ' upper(odefile) ' requires ' msg]); end % Initialize the output function. if haveoutfun PQ = [Q;P]; feval(outfun,[t tfinal],PQ(outputs),'init'); end % Initialize the IRK method method = gniget(options,'Method','G12'); [ns,nm,C,B,BC,AA,E,SM,AM] = coeff_irk2(method); uround = 2.221e-16; B = B * h; BC = BC * h^2; C = C * h; AA = AA * h^2; E(1:ns,1:ns) = E(1:ns,1:ns) * h^2; E(1:ns,(ns+1):(ns+nm)) = E(1:ns,(ns+1):(ns+nm)) * h^2; % ns x ns+nm AM((ns+1):(ns+nm)) = AM((ns+1):(ns+nm)) * h; % ZQ = P * C' + 0.5 * FS * C'.^2; % neq x ns PS = P; EQ = zeros(neq,1); EP = zeros(neq,1); F = zeros(neq,ns); % THE MAIN LOOP done = false; outpoint = 0; addpoint = false; steppos = 0; %while ~done nsteps = round((tfinal - t0)/h); for iter = 1:nsteps if locateevents QOld = Q; POld = P; VOld = V; end tnew = t+h; % if (tnew >= tfinal) % done = true; % end addpoint = false; outpoint = outpoint + 1; if (outpoint >= outstep) outpoint = 0; addpoint = true; end % Start of actual step ynew = y + h*... %% CALL STARTB if (nstps > 0) YH = ZQ * AM(1:ns) + AM(ns+1)*PS + AM(ns+2)*P + Q; ZQ = F * E(:,1:ns)' + FS * E(:,ns+1)'; FS = feval(odefile,t+h,Q,args{:}); F(:,1) = feval(odefile,t+h*SM(nm),YH,args{:}); PS = P; ZQ = ZQ + FS * E(:,ns+2)' + F(:,1) * E(:,ns+nm)' + P * C'; end %% End of STARTB %% Solve nonlinear system dynold = 0; dyno = 1; niter = 0; while (dyno > uround) %% CALL RKNITE QQ = Q * ones(1,ns) + ZQ; if (canvector) F = feval(odefile,t*ones(1,ns) + C',QQ,args{:}); else for i=1:ns F(:,i) = feval(odefile,t + C(i),QQ(:,i),args{:}); end end dyno = 0; NewZQ = P * C' + F * AA'; dyno = sqrt(sum(sum(((ZQ - NewZQ)./max(0.1,abs(Q)*ones(1,ns))).^2)) / (ns*neq)); ZQ = NewZQ; %% END RKNITE niter = niter + 1; if (dynold < dyno) & (dyno < 10*uround) break; end if (niter > maxiter) if (dyno < 1e5*uround) msg = sprintf('Convergence of Gauss-Newton failed after %d iterations.\nObtained error = %0.5g, continuing anyway...',maxiter,dyno); warning(msg); break; else msg = sprintf('Convergence of Gauss-Newton failed after %d iterations.\nObtained error = %0.5g, stopping here.',maxiter,dyno); error(msg); end end dynold = dyno; end %% System solved % LOOP FOR ADVANCING ONE STEP. nstps = nstps + 1; % stats steppos = steppos+1; t = t0+steppos*h; AY = Q; EQ = EQ + h*P + F*BC; Q = AY + EQ; EQ = EQ + (AY - Q); AY = P; EP = EP + F*B; P = AY + EP; EP = EP + (AY - P); if nargout > 0 if outflag == 1 & addpoint nout = nout + 1; if nout > length(tout) tout = [tout; zeros(chunk,1)]; pout = [pout; zeros(chunk,neq)]; qout = [qout; zeros(chunk,neq)]; end tout(nout) = t; qout(nout,:) = Q.'; pout(nout,:) = P.'; end if haveoutfun & addpoint PQ = [Q;P]; if feval(outfun,t,PQ(outputs),'') == 1 return; end end elseif haveoutfun & addpoint PQ = [Q;P]; if feval(outfun,t,PQ(outputs),'') == 1 return; end end if locateevents V = feval(odefile,t,[Q;P],'events',varargin{:}); [nteo,nieo,nqeo,npeo,stop] = gnievents(odefile,VOld,V,POld,P,QOld,Q,[t-h,t],direction,terminal,varargin{:}); teout = [teout;nteo]; ieout = [ieout;nieo]; qeout = [qeout;nqeo]; peout = [peout;npeo]; if stop term = true; break; end end end if (nargout > 0 & outflag == 0) nout = nout + 1; if nout > length(tout) tout = [tout; zeros(1,1)]; pout = [pout; zeros(1,neq)]; qout = [qout; zeros(1,neq)]; end tout(nout) = t; qout(nout,:) = Q.'; pout(nout,:) = P.'; end if haveoutfun feval(outfun,[],[],'done'); end if nargout > 0 tout = tout(1:nout); pout = pout(1:nout,:); qout = qout(1:nout,:); end gnimatlab/gni_lmm2.m100644 12122 12120 21003 7553206145 13453 0ustar hairermathfunction [tout,qout,pout,teout,qeout,peout,ieout] = gni_lmm2(odefile,tspan,y0,options,varargin) %GNI_LMM2 Solves second-order differential equations with a multistep method. % The calling sequence is identical to GNI_IRK2, see HELP GNI_IRK2 for details. % % Details are to be found in the article "GniCodes - Matlab Programs for Geometric % Numerical Integration", available from http://www.unige.ch/math/folks/hairer/ % under item software % % Copyright by Ernst Hairer and Martin Hairer, 10.10.2002. true = logical(1); false = ~true; if nargin == 0 error('Not enough input arguments. See GNI_LMM2.'); elseif ~isstr(odefile) & ~isa(odefile, 'inline') error('First argument must be a single-quoted string. See GNI_LMM2.'); end if nargin == 1 tspan = []; y0 = []; options = []; elseif nargin == 2 y0 = []; options = []; elseif nargin == 3 options = []; elseif ~isempty(options) & ~isa(options,'struct') if (length(tspan) == 1) & (length(y0) == 1) & (min(size(options)) == 1) msg = sprintf('Use gni_lmm2(''%s'',tspan,y0,...) instead.',odefile); error(['Obsolete syntax. ' msg]); else error('Correct syntax is gni_lmm2(''odefile'',tspan,y0,options).'); end end if isstr(odefile) & (exist(odefile) == 4) % a SIMULINK model error('Use the SIM command directly for SIMULINK models.'); return; end % Get default tspan and y0 from odefile if none are specified. if isempty(tspan) | isempty(y0) if (nargout(odefile) < 3) & (nargout(odefile) ~= -1) msg = sprintf('Use gni_lmm2(''%s'',tspan,y0,...) instead.',odefile); error(['No default parameters in ' upper(odefile) '. ' msg]); end [def_tspan,def_y0,def_options] = feval(odefile,[],[],'init',varargin{:}); if isempty(tspan) tspan = def_tspan; end if isempty(y0) y0 = def_y0; end if isempty(options) options = def_options; else options = gniset(def_options,options); end end % Test that tspan is internally consistent. if size(y0,2) ~= 1 y0 = y0'; end ntspan = length(tspan); if ntspan > 2 error('The timespan must consist of one or two numbers.'); end if ntspan == 1 tspan = [0;tspan]; t0 = 0; else t0 = tspan(1); end tfinal = tspan(2); if t0 == tfinal error('The final time must be different from the starting time.'); end tdir = sign(tfinal - t0); outstep = gniget(options,'OutputSteps',1); if (size(outstep) ~= [1 1]) error('The option ''OutputSteps'' must contain a single number'); end outflag = 0; if (outstep > 0) outflag = 1; else outstep = 1; end t = t0; if (mod(length(y0),2) ~= 0) error('The initial value must have 2*N components, where N is the dimension of the system.'); end neq = length(y0)/2; Q = y0(1:neq); P = y0((neq+1):(2*neq)); % By default, hmax is 1/10 of the interval. canvector = strcmp(gniget(options,'Vectorized','off'),'on'); locateevents = strcmp(gniget(options,'Events','off'),'on'); teout = []; qeout = []; peout = []; ieout = []; if locateevents [V,terminal,direction] = feval(odefile,t,[Q;P],'events',varargin{:}); end method = gniget(options,'Method','803'); if nargout > 0 outfun = gniget(options,'OutputFcn'); else outfun = gniget(options,'OutputFcn','odeplot'); end if isempty(outfun) haveoutfun = false; else haveoutfun = true; outputs = gniget(options,'OutputSel',1:(2*neq)); end h = gniget(options,'StepSize'); if (isempty(h)) nsteps = gniget(options,'NumSteps'); if (isempty(nsteps)) h = 1e-2; fprintf('Warning: No initial step size provided, using h=1e-2 instead.\n'); else h = (tfinal - t0)/nsteps; end end if ((~isa(h,'double')) | (size(h) ~= [1 1])) error('The option ''StepSize'' must contain a single number'); end % Allocate memory if we're generating output. chunk = 1; if nargout > 0 chunk = max(ceil(128 / neq),1); tout = zeros(chunk,1); pout = zeros(chunk,neq); qout = zeros(chunk,neq); nout = 1; tout(nout) = t; qout(nout,:) = Q.'; pout(nout,:) = P.'; end if isempty(varargin) args = {}; else args = [{[]} varargin]; end FS = feval(odefile,t,Q,args{:}); [m,n] = size(FS); if n > 1 error([upper(odefile) ' must return a column vector.']) elseif m ~= neq msg = sprintf('an initial condition vector of length 2*%d.',m); error(['Solving ' upper(odefile) ' requires ' msg]); end % Initialize the output function. if haveoutfun QP = [Q;P]; feval(outfun,[t tfinal],QP(outputs),'init'); end % Initialize the IRK method QE = zeros(neq,10); KK = 8; [A,B,C] = coeff_lmm2(method); B=B*h; C=C/h; tdend = t0 + (KK-0.99)*h; stopit = false; if (tdend >= tfinal) stopit = true; tdend = tfinal; end myoptions = gniset(options,'OutputFcn',[],'StepSize',h,... 'OutputSteps',1,'Method','G8'); outpoint = 0; if locateevents [GT,QE,PE,nteout,nqeout,npeout,nieout,term] = gni_irk2(odefile,[t0 tdend],y0,myoptions,varargin{:}); indices = find(nteout < t0+KK/2*h); teout = nteout(indices); qeout = nqeout(indices,:); peout = npeout(indices,:); ieout = nieout(indices); if term stopit = true; end else [GT,QE,PE] = gni_irk2(odefile,[t0 tdend],y0,myoptions,varargin{:}); end if stopit tout = GT; qout = QE; pout = PE; if haveoutfun QP = [QE';PE']; % feval(outfun,GT,[PE(:,outputs)';QE(:,outputs)'],''); feval(outfun,GT,QP(outputs,:),''); feval(outfun,[],[],'done'); end return; end % Generate output for the single-step method. outpoint = 0; t = t0; KK2 = KK/2; for position=2:KK2 addpoint = false; outpoint = outpoint + 1; t = t+h; if (outpoint >= outstep) outpoint = 0; addpoint = true; end if nargout > 0 if outflag == 1 & addpoint nout = nout + 1; if nout > length(tout) tout = [tout; zeros(chunk,1)]; qout = [qout; zeros(chunk,neq)]; pout = [pout; zeros(chunk,neq)]; end tout(nout) = t; qout(nout,:) = QE(position,:); pout(nout,:) = PE(position,:); end if haveoutfun & addpoint QP = [QE(position,:)';PE(position,:)']; if feval(outfun,t,QP(outputs),'') == 1 feval(outfun,[],[],'done'); return; end end elseif haveoutfun & addpoint QP = [QE(position,:)';PE(position,:)']; if feval(outfun,t,QP(outputs),'') == 1 feval(outfun,[],[],'done'); return; end end end if locateevents Q = QE(KK2+1,:)'; P = PE(KK2+1,:)'; V = feval(odefile,t0+KK2*h,[Q;P],'events',varargin{:}); end QE = QE'; Z = (QE(:,2:KK) - QE(:,1:(KK-1)))/h; Z = [Z,zeros(neq,1)]; QE = [QE,zeros(neq,1)]; F = zeros(neq,KK-1); if (canvector) F = feval(odefile,t0+(1:(KK-1))*h,QE(:,2:KK),args{:}); else for i=1:(KK-1) F(:,i) = feval(odefile,t0+i*h,QE(:,i+1),args{:}); end end steppos = KK; t = t0 + steppos*h; U = (Z(:,1:KK2) + Z(:,(KK-1):-1:KK2))*A; % THE MAIN LOOP done = false; addpoint = false; nsteps = round((tfinal - t0)/h)-KK2+1; for iter = 1:nsteps if locateevents QOld = Q; POld = P; VOld = V; end tnew = t+h; addpoint = false; outpoint = outpoint + 1; if (outpoint >= outstep) outpoint = 0; addpoint = true; end % Start of actual step ynew = y + h*... F(:,KK-1) = feval(odefile,t,QE(:,KK),args{:}); U = U + (F(:,1:KK2) + F(:,(KK-1):-1:KK2))*B; Z(:,KK) = zeros(neq,1); SUM = (Z(:,2:(KK2+1)) + Z(:,KK:-1:(KK2+1)))*A; Z(:,KK) = U - SUM; QE(:,KK+1) = QE(:,KK)+h*Z(:,KK); Q = QE(:,KK2+1); P = (QE(:,(KK2+2):(KK+1)) - QE(:,KK2:-1:1))*C; F(:,1:(KK-2)) = F(:,2:(KK-1)); Z(:,2:(KK-1)) = Z(:,3:KK); QE(:,1:KK) = QE(:,2:(KK+1)); % LOOP FOR ADVANCING ONE STEP. if nargout > 0 & outflag == 1 if addpoint nout = nout + 1; if nout > length(tout) tout = [tout; zeros(chunk,1)]; qout = [qout; zeros(chunk,neq)]; pout = [pout; zeros(chunk,neq)]; end tout(nout) = t-KK2*h; qout(nout,:) = Q.'; pout(nout,:) = P.'; end if haveoutfun & addpoint QP = [Q;P]; if feval(outfun,t-KK2*h,QP(outputs),'') == 1 done = true; end end elseif haveoutfun & addpoint QP = [Q;P]; if feval(outfun,t-KK2*h,QP(outputs),'') == 1 done = true; end end if locateevents V = feval(odefile,t-KK2*h,[Q;P],'events',varargin{:}); [nteo,nieo,nqeo,npeo,stop] = gnievents(odefile,VOld,V,POld,P,QOld,Q,[t-(KK2+1)*h,t-KK2*h],direction,terminal,varargin{:}); teout = [teout;nteo]; ieout = [ieout;nieo]; qeout = [qeout;nqeo]; peout = [peout;npeo]; if stop break; end end steppos = steppos+1; t = t0+steppos*h; end if (nargout > 0 & outflag == 0) nout = nout + 1; if nout > length(tout) tout = [tout; zeros(1,1)]; pout = [pout; zeros(1,neq)]; qout = [qout; zeros(1,neq)]; end tout(nout) = t; qout(nout,:) = Q.'; pout(nout,:) = P.'; end if haveoutfun feval(outfun,[],[],'done'); end if nargout > 0 tout = tout(1:nout); qout = qout(1:nout,:); pout = pout(1:nout,:); end gnimatlab/gnievents.m100644 12122 12120 3477 7551505750 13752 0ustar hairermathfunction [teout,ieout,qeout,peout,stop] = gnievents(odefile,E0,E1,P0,P1,Q0,Q1,T,direction,terminal,varargin) H=T(2) - T(1); H2=H*H; DIFF=(Q1-Q0)/H; P0DIF=P0-DIFF; P1DIF=P1-DIFF; teout = []; ieout = []; qeout = []; peout = []; stop = logical(0); LeftSave = []; RightSave = []; TSave = []; LeftCur = E0; RightCur = E1; TCur = T; StateCur = ((E0.*E1) <= 0) & (E1 ~= 0) & (direction .* (E1-E0) >= 0); Level = 0; while (max(StateCur) == 1) Time = (TCur(2)+TCur(1))/2; XM0=Time-T(1); XM1=Time-T(2); POLQ=Q0+XM0*(DIFF+XM1*(P1DIF*XM0+P0DIF*XM1)/H2); POLP=DIFF+(XM0*(2*XM1+XM0)*P1DIF+XM1*(2*XM0+XM1)*P0DIF)/H2; CurVal = feval(odefile,Time,[POLQ;POLP],'events',varargin{:}); StateR = (CurVal.*RightCur) <= 0; StateL = (CurVal.*LeftCur) <= 0; stopit = logical((TCur(1) == Time) | (TCur(2) == Time)); if (StateR == StateCur) LeftCur = CurVal; TCur(1) = Time; elseif (StateL == StateCur) RightCur = CurVal; TCur(2) = Time; else Level = Level+1; LeftSave = [LeftSave,CurVal]; RightSave = [RightSave,RightCur]; TSave = [TSave,[Time;TCur(2)]]; RightCur = CurVal; TCur(2) = Time; StateCur = StateL; end if stopit indices = (find(StateCur)); teout = [teout;Time*ones(size(indices))']; ieout = [ieout;indices]; qeout = [qeout;ones(size(indices))'*POLQ']; peout = [peout;ones(size(indices))'*POLP']; if (max(terminal(indices)) == 1) stop = logical(1); break; end if (Level > 0) TCur = TSave(:,Level); LeftCur = LeftSave(:,Level); RightCur = RightSave(:,Level); StateCur = ((LeftCur.*RightCur) <= 0) & (direction .* (RightCur - LeftCur) >= 0); Level = Level-1; TSave = TSave(:,1:Level); LeftSave = LeftSave(:,1:Level); RightSave = RightSave(:,1:Level); else break; end end end gnimatlab/gniget.m100644 12122 12120 3064 7551506001 13203 0ustar hairermathfunction o = gniget(options,name,default) %GNIGET Get GNI OPTIONS parameters. % VAL = GNIGET(OPTIONS,'NAME') extracts the value of the named property % from the options structure OPTIONS. The syntax is the same as for ODEGET. % % See also GNISET, GNI_IRK2, GNI_COMP, GNI_LMM2. if nargin < 2 error('Not enough input arguments.'); end if ~isempty(options) & ~isa(options,'struct') error('First argument must be an options structure created with GNISET.'); end if isempty(options) if nargin == 3 o = default; else o = []; end return; end Names = [ 'OutputFcn ' 'OutputSel ' 'OutputSteps ' 'Events ' 'Stats ' 'Vectorized ' 'Jacobian ' 'StepSize ' 'NumSteps ' 'MaxIter ' 'PartialFlow ' 'Method ' ]; [m,n] = size(Names); names = lower(Names); lowName = lower(name); j = strmatch(lowName,names); if isempty(j) % if no matches error(sprintf(['Unrecognized property name ''%s''. ' ... 'See GNISET for possibilities.'], name)); elseif length(j) > 1 % if more than one match % Check for any exact matches (in case any names are subsets of others) k = strmatch(lowName,names,'exact'); if length(k) == 1 j = k; else msg = sprintf('Ambiguous property name ''%s'' ', name); msg = [msg '(' deblank(Names(j(1),:))]; for k = j(2:length(j))' msg = [msg ', ' deblank(Names(k,:))]; end msg = sprintf('%s).', msg); error(msg); end end eval(['o = options.' Names(j,:) ';']); if (nargin == 3) & isempty(o) o = default; end gnimatlab/gniset.m100644 12122 12120 12472 7553207007 13250 0ustar hairermathfunction options = gniset(varargin) %GNISET Create/alter GNI OPTIONS structure. % The syntax for GNISET is the same as for ODESET, but the list of possible % properties is different. % %GNISET PROPERTIES % %OutputFcn - Name of installable output function [ string ] % This output function is called by the solver after each time step. When % a solver is called with no output arguments, OutputFcn defaults to % 'odeplot'. Otherwise, OutputFcn defaults to ''. % %OutputSel - Output selection indices [ vector of integers ] % This vector of indices specifies which components of the solution vector % are passed to the OutputFcn. OutputSel defaults to all components. % %OutputSteps - Which steps to output [ integer | 1 ] % This value tells which computed solution points are sent to the output. If % OutputSteps = 10, every 10th solution point is sent to the output. If % OutputSteps = -1, only the first and the last values are sent to the output. % %Vectorized - Vectorized ODE file [ on | {off} ] % Set this property 'on' if the ODE file is coded so that F(t,[y1 y2 ...]) % returns [F(t,y1) F(t,y2) ...]. % %Events - Locate events [ on | {off} ] % Set this property 'on' if the ODE file is coded so that F(t,y,'events') % returns the values of an event function. See ODEFILE. % %StepSize - Step size [ positive scalar ] % The solver will use a fixed step size given by StepSize. It may be % slightly altered if the length of the integration interval is not an % integer multiple of the step size. % %NumSteps - Number of steps [ positive integer ] % The solver will use a fixed step size adjusted to make NumSteps steps. % If StepSize is set, it overrides this option. % %MaxIter - Maximal number of iterations [ integer | {50} ] % Maximal number of Gauss-Newton iterations performed at each step for gni_irk2. % %PartialFlow - Name of a flow function [ string | {'stverl'} ] % The flow function is used by the composition method gni_comp. % %Method - Name of the method [ string | {'817'} ] % For gni_irk2, implemented methods are 'G4', 'G8', and 'G12'. % For gni_comp, implemented methods are '43', '45', '67', '69', % '815', '817', '1033', and '21'. % For gni_lmm2, implemented methods are '801', '802', and '803'. % In all cases the first one or two digits denote the order of the method. % % See also GNIGET, GNI_IRK2, GNI_COMP, GNI_LMM2. % Print out possible values of properties. if (nargin == 0) & (nargout == 0) fprintf(' OutputFcn: [ string ]\n'); fprintf(' OutputSel: [ vector of integers ]\n'); fprintf(' OutputSteps: [ integer ]\n'); fprintf(' Vectorized: [ on | {off} ]\n'); fprintf(' Events: [ on | {off} ]\n'); fprintf(' StepSize: [ positive scalar ]\n'); fprintf(' NumSteps: [ positive integer ]\n'); fprintf(' MaxIter: [ positive integer ]\n'); fprintf(' PartialFlow: [ string ]\n'); fprintf(' Method: [ string ]\n'); fprintf('\n'); return; end Names = [ 'OutputFcn ' 'OutputSel ' 'OutputSteps ' 'Events ' 'Stats ' 'Vectorized ' 'Jacobian ' 'StepSize ' 'NumSteps ' 'MaxIter ' 'PartialFlow ' 'Method ' ]; [m,n] = size(Names); names = lower(Names); options = []; i = 1; while i <= nargin arg = varargin{i}; if isstr(arg) % arg is an option name break; end if ~isempty(arg) % [] is a valid options argument if ~isa(arg,'struct') error(sprintf(['Expected argument %d to be a string property name ' ... 'or an options structure\ncreated with GNISET.'], i)); end if isempty(options) options = arg; else for j = 1:m val = getfield(arg,Names(j,:)); if ~isequal(val,[]) % empty strings '' do overwrite options = setfield(options,Names(j,:),val); end end end end i = i + 1; end if isempty(options) for j = 1:m options = setfield(options,Names(j,:),[]); end end % A finite state machine to parse name-value pairs. if rem(nargin-i+1,2) ~= 0 error('Arguments must occur in name-value pairs.'); end expectval = 0; % start expecting a name, not a value while i <= nargin arg = varargin{i}; if ~expectval if ~isstr(arg) error(sprintf('Expected argument %d to be a string property name.', i)); end lowArg = lower(arg); j = strmatch(lowArg,names); if isempty(j) % if no matches error(sprintf('Unrecognized property name ''%s''.', arg)); elseif length(j) > 1 % if more than one match % Check for any exact matches (in case any names are subsets of others) k = strmatch(lowArg,names,'exact'); if length(k) == 1 j = k; else msg = sprintf('Ambiguous property name ''%s'' ', arg); msg = [msg '(' deblank(Names(j(1),:))]; for k = j(2:length(j))' msg = [msg ', ' deblank(Names(k,:))]; end msg = sprintf('%s).', msg); error(msg); end end expectval = 1; % we expect a value next else options = setfield(options,Names(j,:),arg); expectval = 0; end i = i + 1; end if expectval error(sprintf('Expected value for property ''%s''.', arg)); end gnimatlab/henon.m100644 12122 12120 1236 7553206241 13041 0ustar hairermathfunction [out,out2,out3] = henon(t,q,flag) % Example GNI problem file for solving the Henon-Heiles problem. % Use it as follows: % % [T,Q,P,TE,QE,PE] = gni_meth('henon'); % plot3(QE(:,2),PE(:,1),PE(:,2),'o'); % % With gni_meth replaced by either gni_irk2, gni_lmm2, or gni_comp. if (nargin < 3) | isempty(flag) out(1,:)=-q(1,:).*(1+2*q(2,:)); out(2,:)=-q(2,:).*(1-q(2,:)) - q(1,:).^2; else switch flag case 'init', out = [0 1000]; out2 = [0.18 0.18 0.18 0.18]; out3 = gniset('StepSize',0.2,'Vectorized','on',... 'Events','on','OutputSteps',-1,'OutputFcn','phaseplot'); case 'events', out = [q(1)]; out2 = [0]; out3 = [0]; end end gnimatlab/kepler.m100644 12122 12120 1454 7550620240 13212 0ustar hairermathfunction [out,out2,out3] = kepler(t,q,flag,varargin) % Example GNI problem file for solving the 2D Kepler problem. % Use it as follows: % % gni_meth('kepler',[],[],[],0.7); % % With gni_meth replaced by either gni_irk2, gni_lmm2, or gni_comp % to draw the solution with eccentricity 0.7. if (nargin < 3) | isempty(flag) rad=q(1,:).*q(1,:)+q(2,:).*q(2,:); rad=rad.*sqrt(rad); out(1,:)=-q(1,:)./rad; out(2,:)=-q(2,:)./rad; else switch flag case 'init', if (nargin < 4) ecc = 0.5; else ecc = varargin{1}; end if (ecc < 0) | (ecc >= 1) error('The eccentricity must lie between 0 and 1'); end out = [0 2*pi]; out2 = [1-ecc,0,0,sqrt((1+ecc)/(1-ecc))]; out3 = gniset('NumSteps',50,'Vectorized','on','Events','off',... 'OutputSteps',-1,'OutputFcn','phaseplot'); end end gnimatlab/phaseplot.m100644 12122 12120 5341 7551506004 13730 0ustar hairermathfunction status = phaseplot(t,y,flag) % Example of custom output function. It has to be used in % a GNI options structure as % % gniset('OutputFcn','phaseplot'); % % If 2*N components are given, it plots the N last components % as functions of the N first components. status = 0; chunk = 128; if nargin < 3 | isempty(flag) ud = get(gcf,'UserData'); nt = length(t); chunk = max(chunk,nt); [cols,rows] = size(ud.x); if ud.i + nt > rows ud.x = [ud.x, zeros(ud.N,chunk)]; ud.y = [ud.y, zeros(ud.N,chunk)]; end ud.x(1:ud.N,(ud.i+1):(ud.i+nt)) = y(1:ud.N,:); ud.y(1:ud.N,(ud.i+1):(ud.i+nt)) = y((ud.N+1):(2*ud.N),:); ud.i = ud.i + nt; set(gcf,'UserData',ud); if ud.stop == 1 status = 1; phaseplot([],[],'done'); else xlim = get(gca,'xlim'); ylim = get(gca,'ylim'); if (ud.i == nt+1) | (min(y(1:ud.N,1:nt)) < xlim(1)) | ... (xlim(2) < max(y(1:ud.N,1:nt))) | ... (min(y((ud.N+1):(2*ud.N),1:nt)) < ylim(1)) | ... (ylim(2) < max(y((ud.N+1):(2*ud.N),1:nt))) for i=1:ud.N set(ud.lines(i),'Xdata',ud.x(i,1:ud.i),'Ydata',ud.y(i,1:ud.i)); end else for i=1:ud.N set(ud.line(i),'Xdata',ud.x(i,1:ud.i),'Ydata',ud.y(i,1:ud.i)); end end end else switch(flag) case 'init' ud = []; ud.N = floor(size(y) / 2); ud.N = ud.N(1); ud.x = zeros(ud.N,chunk); ud.y = zeros(ud.N,chunk); ud.i = 1; for i=1:ud.N ud.x(i,1) = y(i); ud.y(i,1) = y(i+ud.N); end f = figure(gcf); hold on colr = ['b','r','g','m','c']; for i=1:ud.N ud.line(i) = plot(ud.x(i,1),ud.y(i,1),... strcat(colr(mod(i-1,5)+1),'o-'),'EraseMode','none'); ud.lines(i) = plot(ud.x(i,1),ud.y(i,1),... strcat(colr(mod(i-1,5)+1),'o-')); end set(gca,'XLimMode','auto'); set(gca,'YLimMode','auto'); hold off % The STOP button. h = findobj(f,'Tag','stop'); if isempty(h) ud.stop = 0; pos = get(0,'DefaultUicontrolPosition'); pos(1) = pos(1) - 15; pos(2) = pos(2) - 15; str = 'ud=get(gcf,''UserData''); ud.stop=1; set(gcf,''UserData'',ud);'; uicontrol( ... 'Style','push', ... 'String','Stop', ... 'Position',pos, ... 'Callback',str, ... 'Tag','stop'); else set(h,'Visible','on'); if ishold oud = get(f,'UserData'); ud.stop = oud.stop; else ud.stop = 0; end end set(f,'UserData',ud); case 'done' f = gcf; ud = get(f,'UserData'); ud.x = ud.x(1:ud.N,1:ud.i); ud.y = ud.y(1:ud.N,1:ud.i); set(f,'UserData',ud); for i=1:ud.N set(ud.line(i),'Xdata',ud.x(i,1:ud.i),'Ydata',ud.y(i,1:ud.i)); end set(findobj(f,'Tag','stop'),'Visible','off'); refresh; end end drawnow; gnimatlab/rattwo.m100755 12122 12120 1307 7551505770 13263 0ustar hairermathfunction [outP,outQ] = rattwo(t,P,Q,ode,ha,hb,hc,first,last,flags,varargin) if isempty(flags) F = feval(ode,t,Q,varargin{:}); EP = P - ha*F; EQ = Q + hb*EP; EE1 = EQ(1:3)'*EQ(1:3); EQ1 = EQ(1:3)'*Q(1:3); EE2 = EQ(4:6)'*EQ(4:6); EQ2 = EQ(4:6)'*Q(4:6); BET1 = 1 - EE1; ALAM1 = -BET1/(hb*(EQ1+sqrt(BET1+EQ1^2))); BET2 = 1 - EE2; ALAM2 = -BET2/(hb*(EQ2+sqrt(BET2+EQ2^2))); outP = EP - [ALAM1*Q(1:3);ALAM2*Q(4:6)]; outQ = Q + hb*outP; if (last) F = feval(ode,t,outQ,varargin{:}); outP = outP - hc*F; AMU1 = sum(outP(1:3).*outQ(1:3)); AMU2 = sum(outP(4:6).*outQ(4:6)); outP = outP - [AMU1*outQ(1:3);AMU2*outQ(4:6)]; end return; else switch flags case 'init', case 'done', end end gnimatlab/sphereplot.m100644 12122 12120 5224 7551505773 14132 0ustar hairermathfunction status = sphereplot(t,y,flag) % Example of custom output function. It has to be used in % a GNI options structure as % % gniset('OutputFcn','sphereplot'); % % It requires to have 3*N components given. It interprets them % as N cartesian coordinate triplets that are supposed to lie % on the unit sphere. status = 0; chunk = 128; if nargin < 3 | isempty(flag) ud = get(gcf,'UserData'); nt = length(t); chunk = max(chunk,nt); [cols,rows] = size(ud.x); if ud.i + 1 > rows ud.x = [ud.x, zeros(ud.N,chunk)]; ud.y = [ud.y, zeros(ud.N,chunk)]; ud.z = [ud.z, zeros(ud.N,chunk)]; end ud.i = ud.i + 1; ud.x(1:ud.N,ud.i) = y(1:3:(3*ud.N)); ud.y(1:ud.N,ud.i) = y(2:3:(3*ud.N)); ud.z(1:ud.N,ud.i) = y(3:3:(3*ud.N)); set(gcf,'UserData',ud); if ud.stop == 1 status = 1; sphereplot([],[],'done'); else for i=1:ud.N set(ud.line(i),'Xdata',ud.x(i,1:ud.i),'Ydata',ud.y(i,1:ud.i),... 'Zdata',ud.z(i,1:ud.i)); end end else switch(flag) case 'init' ud = []; ud.N = floor(size(y) / 3); ud.N = ud.N(1); ud.x = zeros(ud.N,chunk); ud.y = zeros(ud.N,chunk); ud.z = zeros(ud.N,chunk); ud.i = 1; for i=1:ud.N ud.x(i,1) = y(3*i-2); ud.y(i,1) = y(3*i-1); ud.z(i,1) = y(3*i); end f = figure(gcf); colormap([0.8 0.8 1]); [X,Y,Z] = sphere(100); sph = surfl(X,Y,Z,'light'); set(sph(1),'EdgeColor','none'); hold on colr = ['b','r','g','m','c']; for i=1:ud.N ud.line(i) = plot3(ud.x(i,1),ud.y(i,1),ud.z(i,1),... strcat(colr(mod(i-1,5)+1),'o-'),'EraseMode','none'); end hold off set(gca,'XLim',[-1 1]); set(gca,'YLim',[-1 1]); set(gca,'ZLim',[-1 1]); % The STOP button. h = findobj(f,'Tag','stop'); if isempty(h) ud.stop = 0; pos = get(0,'DefaultUicontrolPosition'); pos(1) = pos(1) - 15; pos(2) = pos(2) - 15; str = 'ud=get(gcf,''UserData''); ud.stop=1; set(gcf,''UserData'',ud);'; uicontrol( ... 'Style','push', ... 'String','Stop', ... 'Position',pos, ... 'Callback',str, ... 'Tag','stop'); else set(h,'Visible','on'); if ishold oud = get(f,'UserData'); ud.stop = oud.stop; else ud.stop = 0; end end set(f,'UserData',ud); case 'done' f = gcf; ud = get(f,'UserData'); ud.x = ud.x(1:ud.N,1:ud.i); ud.y = ud.y(1:ud.N,1:ud.i); ud.z = ud.z(1:ud.N,1:ud.i); set(f,'UserData',ud); for i=1:ud.N set(ud.line(i),'Xdata',ud.x(i,1:ud.i),'Ydata',ud.y(i,1:ud.i)... ,'Zdata',ud.z(i,1:ud.i)); end set(findobj(f,'Tag','stop'),'Visible','off'); refresh; end end drawnow; gnimatlab/stverl.m100755 12122 12120 672 7551505756 13252 0ustar hairermathfunction [outP,outQ] = stverl(t,P,Q,ode,ha,hb,hc,first,last,flags,varargin) global EP; global EQ; if isempty(flags) if (first) EQ = EQ + ha*P; outQ = Q + EQ; EQ = EQ + (Q - outQ); Q = outQ; end F = feval(ode,t,Q,varargin{:}); EP = EP + hb*F; outP = P + EP; EP = EP + (P - outP); EQ = EQ + hc*outP; outQ = Q + EQ; EQ = EQ + (Q - outQ); return; else switch flags case 'init', EQ = 0; EP = 0; case 'done', end end gnimatlab/twobodysphere.m100644 12122 12120 2134 7552754654 14645 0ustar hairermathfunction [out,out2,out3] = twobodysphere(t,q,flag) % Example GNI problem file for solving the two-body problem % constrained to the sphere of radius one. % Use it as follows: % % gni_comp('twobodysphere'); % % Only gni_comp works, since this problem uses a special basic % method and is therefore not written in the standard way. This problem % also uses the custom output function SPHEREPLOT. if (nargin < 3) | isempty(flag) prod = q(1:3)'*q(4:6); out = -q([4:6,1:3])/(1-prod^2)^(3/2); else switch flag case 'init', out = [0 10]; phi = [1.3 -2.1]; theta = [2.1 -1.1]; out2([1 4]) = cos(phi).*sin(theta); out2([2 5]) = sin(phi).*sin(theta); out2([3 6]) = cos(theta); dphi = [1.2 0.1]; dtheta = [0.1 -0.5]; out2([7 10]) = -dphi.*sin(phi).*sin(theta) + dtheta.*cos(phi).*cos(theta); out2([8 11]) = dphi.*cos(phi).*sin(theta) + dtheta.*sin(phi).*cos(theta); out2([9 12]) = -dtheta.*sin(theta); out3 = gniset('StepSize',0.02,'Vectorized','off','Events','off',... 'PartialFlow','rattwo','OutputFcn','sphereplot',... 'OutputSteps',5,'OutputSel',[1:6]); end end