source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int/rosenbrock.m @ 3997

Last change on this file since 3997 was 2696, checked in by kanani, 6 years ago

Merge of branch palm4u into trunk

File size: 31.2 KB
Line 
1function [T, Y, RCNTRL, ICNTRL, RSTAT, ISTAT] = ...
2    Rosenbrock(Function, Tspan, Y0, Options, RCNTRL, ICNTRL)
3%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
4%  Rosenbrock - Implementation of several Rosenbrock methods:
5%               * Ros2                                               
6%               * Ros3                                               
7%               * Ros4                                               
8%               * Rodas3                                             
9%               * Rodas4                                             
10%                                                                     
11%    Solves the system y'=F(t,y) using a Rosenbrock method defined by:
12%                                                                     
13%     G = 1/(H*gamma(1)) - Jac(t0,Y0)                                 
14%     T_i = t0 + Alpha(i)*H                                           
15%     Y_i = Y0 + \sum_{j=1}^{i-1} A(i,j)*K_j
16%     G * K_i = Fun( T_i, Y_i ) + \sum_{j=1}^S C(i,j)/H * K_j +
17%         gamma(i)*dF/dT(t0, Y0)
18%     Y1 = Y0 + \sum_{j=1}^S M(j)*K_j
19%
20%    For details on Rosenbrock methods and their implementation consult:
21%      E. Hairer and G. Wanner
22%      "Solving ODEs II. Stiff and differential-algebraic problems".
23%      Springer series in computational mathematics, Springer-Verlag, 1996.
24%    The codes contained in the book inspired this implementation.
25%
26%    MATLAB implementation (C) John C. Linford (jlinford@vt.edu).     
27%    Virginia Polytechnic Institute and State University             
28%    November, 2009   
29%
30%    Based on the Fortran90 implementation (C) Adrian Sandu, August 2004
31%    and revised by Philipp Miehe and Adrian Sandu, May 2006.
32%    Virginia Polytechnic Institute and State University
33%    Contact: sandu@cs.vt.edu
34%
35%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
36%  Input Arguments :
37%  The first four arguments are similar to the input arguments of
38%  MATLAB's ODE solvers
39%      Function  - A function handle for the ODE function
40%      Tspan     - The time space to integrate
41%      Y0        - Initial value
42%      Options   - ODE solver options created by odeset():           
43%                  AbsTol, InitialStep, Jacobian, MaxStep, and RelTol
44%                  are considered.  Other options are ignored.       
45%                  'Jacobian' must be set.                           
46%      RCNTRL    - real value input parameters (explained below)
47%      ICNTRL    - integer input parameters (explained below)
48%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
49%  Output Arguments:
50%  The first two arguments are similar to the output arguments of
51%  MATLAB's ODE solvers
52%      T         - A vector of final integration times.
53%      Y         - A matrix of function values.  Y(T(i),:) is the value of
54%                  the function at the ith output time.
55%      RCNTRL    - real value input parameters (explained below)
56%      ICNTRL    - integer input parameters (explained below)
57%      RSTAT     - real output parameters (explained below)
58%      ISTAT     - integer output parameters (explained below)
59%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
60%
61%    RCNTRL and ICNTRL on input and output:
62%
63%    Note: For input parameters equal to zero the default values of the
64%       corresponding variables are used.
65%
66%    ICNTRL(1) = 1: F = F(y)   Independent of T (AUTONOMOUS)
67%              = 0: F = F(t,y) Depends on T (NON-AUTONOMOUS)
68%
69%    ICNTRL(2) = 0: AbsTol, RelTol are N-dimensional vectors
70%              = 1: AbsTol, RelTol are scalars
71%
72%    ICNTRL(3)  -> selection of a particular Rosenbrock method
73%        = 0 :    Rodas3 (default)
74%        = 1 :    Ros2
75%        = 2 :    Ros3
76%        = 3 :    Ros4
77%        = 4 :    Rodas3
78%        = 5 :    Rodas4
79%
80%    ICNTRL(4)  -> maximum number of integration steps
81%        For ICNTRL(4)=0) the default value of 100000 is used
82%
83%    RCNTRL(1)  -> Hmin, lower bound for the integration step size
84%          It is strongly recommended to keep Hmin = ZERO
85%    RCNTRL(2)  -> Hmax, upper bound for the integration step size
86%    RCNTRL(3)  -> Hstart, starting value for the integration step size
87%
88%    RCNTRL(4)  -> FacMin, lower bound on step decrease factor (default=0.2)
89%    RCNTRL(5)  -> FacMax, upper bound on step increase factor (default=6)
90%    RCNTRL(6)  -> FacRej, step decrease factor after multiple rejections
91%                          (default=0.1)
92%    RCNTRL(7)  -> FacSafe, by which the new step is slightly smaller
93%         than the predicted value  (default=0.9)
94%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
95%
96%    RSTAT and ISTAT on output:
97%
98%    ISTAT(1)  -> No. of function calls
99%    ISTAT(2)  -> No. of jacobian calls
100%    ISTAT(3)  -> No. of steps
101%    ISTAT(4)  -> No. of accepted steps
102%    ISTAT(5)  -> No. of rejected steps (except at very beginning)
103%    ISTAT(6)  -> No. of LU decompositions
104%    ISTAT(7)  -> No. of forward/backward substitutions
105%    ISTAT(8)  -> No. of singular matrix decompositions
106%
107%    RSTAT(1)  -> Texit, the time corresponding to the
108%                     computed Y upon return
109%    RSTAT(2)  -> Hexit, last accepted step before exit
110%    RSTAT(3)  -> Hnew, last predicted step (not yet taken)
111%                   For multiple restarts, use Hnew as Hstart
112%                     in the subsequent run
113%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
114
115%~~~>  Statistics on the work performed by the Rosenbrock method
116global Nfun Njac Nstp Nacc Nrej Ndec Nsol Nsng Ntexit Nhexit Nhnew
117global ISTATUS RSTATUS
118
119% Parse ODE options
120AbsTol = odeget(Options, 'AbsTol');
121if isempty(AbsTol)
122    AbsTol = 1.0e-3;
123end
124Hstart = odeget(Options, 'InitialStep');
125if isempty(Hstart)
126    Hstart = 0;
127end
128Jacobian = odeget(Options, 'Jacobian');
129if isempty(Jacobian)
130    error('A Jacobian function is required.');
131end
132Hmax = odeget(Options, 'MaxStep');
133if isempty(Hmax)
134    Hmax = 0;
135end
136RelTol = odeget(Options, 'RelTol');
137if isempty(RelTol)
138    RelTol = 1.0e-4;
139end
140
141% Initialize statistics
142Nfun=1; Njac=2; Nstp=3; Nacc=4; Nrej=5; Ndec=6; Nsol=7; Nsng=8;
143Ntexit=1; Nhexit=2; Nhnew=3;
144RSTATUS = zeros(20,1);
145ISTATUS = zeros(20,1);
146
147% Get problem size
148steps = length(Tspan);
149N = length(Y0);
150
151% Initialize tolerances
152ATOL = ones(N,1)*AbsTol;
153RTOL = ones(N,1)*RelTol;
154
155% Initialize step
156RCNTRL(2) = max(0, Hmax);
157RCNTRL(3) = max(0, Hstart);
158
159% Integrate
160Y = zeros(steps,N);
161T = zeros(steps,1);
162Y(1,:) = Y0;
163T(1) = Tspan(1);
164for i=2:steps
165    [T(i), Y(i,:), IERR] = ...
166        RosenbrockIntegrator(N, Y(i-1,:), T(i-1), Tspan(i), ...
167            Function, Jacobian, ATOL, RTOL, RCNTRL, ICNTRL);
168
169    if IERR < 0
170        error(['Rosenbrock exited with IERR=',num2str(IERR)]);
171    end
172
173end
174
175ISTAT = ISTATUS;
176RSTAT = RSTATUS;
177
178return
179
180%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
181function [ T, Y, IERR ] = RosenbrockIntegrator(N, Y, Tstart, Tend, ...
182    OdeFunction, OdeJacobian, AbsTol, RelTol, RCNTRL, ICNTRL)
183%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
184% Advances the ODE system defined by OdeFunction, OdeJacobian and Y
185% from time=Tstart to time=Tend
186%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
187
188%~~~>  Statistics on the work performed by the Rosenbrock method
189global Nfun Njac Nstp Nacc Nrej Ndec Nsol Nsng Ntexit Nhexit Nhnew Ntotal
190global ISTATUS RSTATUS
191
192%~~~>   Parameters
193global ZERO ONE DeltaMin
194ZERO = 0.0; ONE = 1.0; DeltaMin = 1.0E-5;
195
196%~~~>  Initialize statistics
197ISTATUS(1:8) = 0;
198RSTATUS(1:3) = ZERO;
199
200%~~~>  Autonomous or time dependent ODE. Default is time dependent.
201Autonomous = ~(ICNTRL(1) == 0);
202
203%~~~>  For Scalar tolerances (ICNTRL(2).NE.0)  the code uses AbsTol(1) and RelTol(1)
204%   For Vector tolerances (ICNTRL(2) == 0) the code uses AbsTol(1:N) and RelTol(1:N)
205if ICNTRL(2) == 0
206    VectorTol = true;
207    UplimTol  = N;
208else
209    VectorTol = false;
210    UplimTol  = 1;
211end
212
213%~~~>   Initialize the particular Rosenbrock method selected
214switch (ICNTRL(3))
215    case (1)
216        ros_Param = Ros2;
217    case (2)
218        ros_Param = Ros3;
219    case (3)
220        ros_Param = Ros4;
221    case {0,4}
222        ros_Param = Rodas3;
223    case (5)
224        ros_Param = Rodas4;
225    otherwise
226        disp(['Unknown Rosenbrock method: ICNTRL(3)=',num2str(ICNTRL(3))]);
227        IERR = ros_ErrorMsg(-2,Tstart,ZERO);
228        return
229end
230ros_S     = ros_Param{1};
231rosMethod = ros_Param{2};
232ros_A     = ros_Param{3};
233ros_C     = ros_Param{4};
234ros_M     = ros_Param{5};
235ros_E     = ros_Param{6};
236ros_Alpha = ros_Param{7};
237ros_Gamma = ros_Param{8};
238ros_ELO   = ros_Param{9};
239ros_NewF  = ros_Param{10};
240ros_Name  = ros_Param{11};
241
242%~~~>   The maximum number of steps admitted
243if ICNTRL(4) == 0
244    Max_no_steps = 200000;
245elseif ICNTRL(4) > 0
246    Max_no_steps=ICNTRL(4);
247else
248    disp(['User-selected max no. of steps: ICNTRL(4)=',num2str(ICNTRL(4))]);
249    IERR = ros_ErrorMsg(-1,Tstart,ZERO);
250    return
251end
252
253%~~~>  Unit roundoff (1+Roundoff>1)
254Roundoff = eps;
255
256%~~~>  Lower bound on the step size: (positive value)
257if RCNTRL(1) == ZERO
258    Hmin = ZERO;
259elseif RCNTRL(1) > ZERO
260    Hmin = RCNTRL(1);
261else
262    disp(['User-selected Hmin: RCNTRL(1)=', num2str(RCNTRL(1))]);
263    IERR = ros_ErrorMsg(-3,Tstart,ZERO);
264    return
265end
266
267%~~~>  Upper bound on the step size: (positive value)
268if RCNTRL(2) == ZERO
269    Hmax = abs(Tend-Tstart);
270elseif RCNTRL(2) > ZERO
271    Hmax = min(abs(RCNTRL(2)),abs(Tend-Tstart));
272else
273    disp(['User-selected Hmax: RCNTRL(2)=', num2str(RCNTRL(2))]);
274    IERR = ros_ErrorMsg(-3,Tstart,ZERO);
275    return
276end
277
278%~~~>  Starting step size: (positive value)
279if RCNTRL(3) == ZERO
280    Hstart = max(Hmin,DeltaMin);
281elseif RCNTRL(3) > ZERO
282    Hstart = min(abs(RCNTRL(3)),abs(Tend-Tstart));
283else
284    disp(['User-selected Hstart: RCNTRL(3)=', num2str(RCNTRL(3))]);
285    IERR = ros_ErrorMsg(-3,Tstart,ZERO);
286    return
287end
288
289%~~~>  Step size can be changed s.t.  FacMin < Hnew/Hold < FacMax
290if RCNTRL(4) == ZERO
291    FacMin = 0.2;
292elseif RCNTRL(4) > ZERO
293    FacMin = RCNTRL(4);
294else
295    disp(['User-selected FacMin: RCNTRL(4)=', num2str(RCNTRL(4))]);
296    IERR = ros_ErrorMsg(-4,Tstart,ZERO);
297    return
298end
299if RCNTRL(5) == ZERO
300    FacMax = 6.0;
301elseif RCNTRL(5) > ZERO
302    FacMax = RCNTRL(5);
303else
304    disp(['User-selected FacMax: RCNTRL(5)=', num2str(RCNTRL(5))]);
305    IERR = ros_ErrorMsg(-4,Tstart,ZERO);
306    return
307end
308
309%~~~>   FacRej: Factor to decrease step after 2 succesive rejections
310if RCNTRL(6) == ZERO
311    FacRej = 0.1;
312elseif RCNTRL(6) > ZERO
313    FacRej = RCNTRL(6);
314else
315    disp(['User-selected FacRej: RCNTRL(6)=', num2str(RCNTRL(6))]);
316    IERR = ros_ErrorMsg(-4,Tstart,ZERO);
317    return
318end
319
320%~~~>   FacSafe: Safety Factor in the computation of new step size
321if RCNTRL(7) == ZERO
322    FacSafe = 0.9;
323elseif RCNTRL(7) > ZERO
324    FacSafe = RCNTRL(7);
325else
326    disp(['User-selected FacSafe: RCNTRL(7)=', num2str(RCNTRL(7))]);
327    IERR = ros_ErrorMsg(-4,Tstart,ZERO);
328    return
329end
330
331%~~~>  Check if tolerances are reasonable
332for i=1:UplimTol
333    if  (AbsTol(i) <= ZERO) || (RelTol(i) <= 10.0*Roundoff) || (RelTol(i) >= 1.0)
334        disp([' AbsTol(',i,') = ',num2str(AbsTol(i))]);
335        disp([' RelTol(',i,') = ',num2str(RelTol(i))]);
336        IERR = ros_ErrorMsg(-5,Tstart,ZERO);
337        return
338    end
339end
340
341%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
342%   Template for the implementation of a generic Rosenbrock method
343%      defined by ros_S (no of stages)
344%      and its coefficients ros_{A,C,M,E,Alpha,Gamma}
345%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
346
347% ~~~~ Local variables
348Ynew  = zeros(N,1);
349Fcn0  = zeros(N,1);
350Fcn   = zeros(N,1);
351K     = zeros(N,ros_S);
352dFdT  = zeros(N,1);
353Jac0  = sparse(N,N);
354Ghimj = sparse(N,N);
355Yerr  = zeros(N,1);
356Pivot = zeros(N,1);
357
358%~~~>  Initial preparations
359T = Tstart;
360RSTATUS(Nhexit) = ZERO;
361H = min(max(abs(Hmin), abs(Hstart)), abs(Hmax));
362if abs(H) <= 10.0*Roundoff
363    H = DeltaMin;
364end
365
366if Tend >= Tstart
367    Direction = +1;
368else
369    Direction = -1;
370end
371H = Direction*H;
372
373RejectLastH = false;
374RejectMoreH = false;
375
376%~~~> Time loop begins below
377
378while ( (Direction > 0) && ((T-Tend)+Roundoff <= ZERO) || ...
379        (Direction < 0) && ((Tend-T)+Roundoff <= ZERO) )
380
381    if ISTATUS(Nstp) > Max_no_steps   % Too many steps
382        IERR = ros_ErrorMsg(-6,T,H);
383        return
384    end
385    if ((T+0.1*H) == T) || (H <= Roundoff)   % Step size too small
386        IERR = ros_ErrorMsg(-7,T,H);
387        return
388    end
389
390    %~~~>  Limit H if necessary to avoid going beyond Tend
391    H = min(H, abs(Tend-T));
392
393    %~~~>   Compute the function at current time
394    Fcn0 = OdeFunction(T, Y);
395    ISTATUS(Nfun) = ISTATUS(Nfun) + 1;
396
397    %~~~>  Compute the function derivative with respect to T
398    if ~Autonomous
399        dFdT = ros_FunTimeDerivative (T, OdeFunction, Roundoff, Y, Fcn0);
400    end
401
402    %~~~>   Compute the Jacobian at current time
403    Jac0 = OdeJacobian(T,Y);
404    ISTATUS(Njac) = ISTATUS(Njac) + 1;
405
406    %~~~>  Repeat step calculation until current step accepted
407    while(true)
408
409        [H, Ghimj, Pivot, Singular] = ...
410            ros_PrepareMatrix(N, H, Direction, ros_Gamma(1), Jac0);
411
412        % Not calculating LU decomposition in ros_PrepareMatrix anymore
413        % so don't need to check if the matrix is singular
414        %
415        %if Singular % More than 5 consecutive failed decompositions
416        %    IERR = ros_ErrorMsg(-8,T,H);
417        %    return
418        %end
419
420        % For the 1st istage the function has been computed previously
421        Fcn = Fcn0;
422        K(:,1) = Fcn;
423        if (~Autonomous) && (ros_Gamma(1) ~= ZERO)
424            HG = Direction*H*ros_Gamma(1);
425            K(:,1) = K(:,1) + HG * dFdT;
426        end
427        K(:,1) = ros_Solve(Ghimj, Pivot, K(:,1));
428       
429        %~~~>   Compute the remaining stages
430        for istage=2:ros_S
431
432            % istage>1 and a new function evaluation is needed
433            if ros_NewF(istage)
434                Ynew = Y;
435                for j=1:istage-1
436                    Ynew = Ynew + ros_A((istage-1)*(istage-2)/2+j)*K(:,j)';
437                end
438                Tau = T + ros_Alpha(istage)*Direction*H;
439                Fcn = OdeFunction(Tau,Ynew);
440                ISTATUS(Nfun) = ISTATUS(Nfun) + 1;
441            end
442
443            K(:,istage) = Fcn;
444            for j=1:istage-1
445                HC = ros_C((istage-1)*(istage-2)/2+j)/(Direction*H);
446                K(:,istage) = K(:,istage) + HC * K(:,j);
447            end
448            if (~Autonomous) && (ros_Gamma(istage) ~= ZERO)
449                HG = Direction*H*ros_Gamma(istage);
450                K(:,istage) = K(:,istage) + HG * dFdT;
451            end
452
453            % Linear system is solved here with MATLAB '\' operator
454            % instead of LU decomposition.
455            K(:,istage) = ros_Solve(Ghimj, Pivot, K(:,istage));
456
457        end
458
459        %~~~>  Compute the new solution
460        Ynew = Y;
461        for j=1:ros_S
462            Ynew = Ynew + ros_M(j) * K(:,j)';
463        end
464
465        %~~~>  Compute the error estimation
466        Yerr = zeros(N,1);
467        for j=1:ros_S
468            Yerr = Yerr + ros_E(j) * K(:,j);
469        end
470        Err = ros_ErrorNorm(N, Y, Ynew, Yerr, AbsTol, RelTol, VectorTol);
471
472        %~~~> New step size is bounded by FacMin <= Hnew/H <= FacMax
473        Fac  = min(FacMax, max(FacMin, FacSafe/(Err^(ONE/ros_ELO))));
474        Hnew = H*Fac;
475
476        %~~~>  Check the error magnitude and adjust step size
477        ISTATUS(Nstp) = ISTATUS(Nstp) + 1;
478        if  (Err <= ONE) || (H <= Hmin)   %~~~> Accept step
479            ISTATUS(Nacc) = ISTATUS(Nacc) + 1;
480            Y = Ynew;
481            T = T + Direction*H;
482            Hnew = max(Hmin, min(Hnew, Hmax));
483            if RejectLastH  % No step size increase after a rejected step
484                Hnew = min(Hnew, H);
485            end
486            RSTATUS(Nhexit) = H;
487            RSTATUS(Nhnew)  = Hnew;
488            RSTATUS(Ntexit) = T;
489            RejectLastH = false;
490            RejectMoreH = false;
491            H = Hnew;
492            break % EXIT THE LOOP: WHILE STEP NOT ACCEPTED
493        else           %~~~> Reject step
494            if RejectMoreH
495                Hnew = H*FacRej;
496            end
497            RejectMoreH = RejectLastH;
498            RejectLastH = true;
499            H = Hnew;
500            if ISTATUS(Nacc) >= 1
501                ISTATUS(Nrej) = ISTATUS(Nrej) + 1;
502            end
503        end % Err <= 1
504
505    end % UntilAccepted
506
507end % TimeLoop
508
509%~~~> Succesful exit
510IERR = 1;  %~~~> The integration was successful
511
512return
513
514
515%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
516function [ ErrNorm ] = ros_ErrorNorm(N, Y, Ynew, Yerr, AbsTol, RelTol, VectorTol)
517%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
518%~~~> Computes the "scaled norm" of the error vector Yerr
519%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
520
521Err = 0;
522for i=1:N
523    Ymax = max(abs(Y(i)), abs(Ynew(i)));
524    if VectorTol
525        Scale = AbsTol(i) + RelTol(i)*Ymax;
526    else
527        Scale = AbsTol(1) + RelTol(1)*Ymax;
528    end
529    Err = Err + (Yerr(i)/Scale)^2;
530end
531Err = sqrt(Err/N);
532ErrNorm = max(Err,1.0e-10);
533
534return
535
536%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
537function [ dFdT ] = ros_FunTimeDerivative (T, OdeFunction, Roundoff, Y, Fcn0)
538%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
539%~~~> The time partial derivative of the function by finite differences
540%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
541
542%~~~>  Statistics on the work performed by the Rosenbrock method
543global Nfun ISTATUS
544
545%~~~>   Parameters
546global ZERO ONE DeltaMin
547
548Delta = sqrt(Roundoff) * max(DeltaMin, abs(T));
549dFdT = OdeFunction(T+Delta,Y);
550ISTATUS(Nfun) = ISTATUS(Nfun) + 1;
551dFdT = (dFdT - Fcn0) / Delta;
552
553return
554
555
556%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
557function [H, Ghimj, Pivot, Singular] = ...
558    ros_PrepareMatrix(N, H, Direction, gam, Jac0)
559%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
560%  Prepares the LHS matrix for stage calculations
561%  1.  Construct Ghimj = 1/(H*ham) - Jac0
562%      "(Gamma H) Inverse Minus Jacobian"
563%  2.  LU decomposition not performed here because
564%      MATLAB solves the linear system with '\'.
565%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
566
567Ghimj = -Jac0;
568ghinv = 1.0/(Direction*H*gam);
569for i=1:N
570    Ghimj(i,i) = Ghimj(i,i)+ghinv;
571end
572
573Pivot(1) = 1;
574Singular = false;
575
576return
577
578
579%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
580function [ IERR ] = ros_ErrorMsg(Code, T, H)
581%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
582%    Handles all error messages
583%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
584
585IERR = Code;
586
587disp('Forced exit from Rosenbrock due to the following error:');
588
589switch (Code)
590    case (-1)
591        disp('--> Improper value for maximal no of steps');
592    case (-2)
593        disp('--> Selected Rosenbrock method not implemented');
594    case (-3)
595        disp('--> Hmin/Hmax/Hstart must be positive');
596    case (-4)
597        disp('--> FacMin/FacMax/FacRej must be positive');
598    case (-5)
599        disp('--> Improper tolerance values');
600    case (-6)
601        disp('--> No of steps exceeds maximum bound');
602    case (-7)
603        disp('--> Step size too small: T + 10*H = T or H < Roundoff');
604    case (-8)
605        disp('--> Matrix is repeatedly singular');
606    otherwise
607        disp(['Unknown Error code: ', num2str(Code)]);
608end
609
610disp(['T=',num2str(T),' and H=',num2str(H)]);
611
612return
613
614
615%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
616function [X] = ros_Solve(JVS, Pivot, X)
617%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
618%  Template for the forward/backward substitution
619%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
620
621global Nsol ISTATUS
622
623X = JVS\X;
624Pivot(1) = 1;
625ISTATUS(Nsol) = ISTATUS(Nsol) + 1;
626
627return
628
629
630%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
631function [ params ] = Ros2()
632%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
633% --- AN L-STABLE METHOD, 2 stages, order 2
634%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
635
636g = 1.0 + 1.0/sqrt(2.0);
637rosMethod = 1;
638%~~~> Name of the method
639ros_Name = 'ROS-2';
640%~~~> Number of stages
641ros_S = 2;
642
643%~~~> The coefficient matrices A and C are strictly lower triangular.
644%   The lower triangular (subdiagonal) elements are stored in row-wise order:
645%   A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
646%   The general mapping formula is:
647%       A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
648%       C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
649
650ros_A(1) = (1.0)/g;
651ros_C(1) = (-2.0)/g;
652%~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
653%   or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
654ros_NewF(1) = true;
655ros_NewF(2) = true;
656%~~~> M_i = Coefficients for new step solution
657ros_M(1)= (3.0)/(2.0*g);
658ros_M(2)= (1.0)/(2.0*g);
659% E_i = Coefficients for error estimator
660ros_E(1) = 1.0/(2.0*g);
661ros_E(2) = 1.0/(2.0*g);
662%~~~> ros_ELO = estimator of local order - the minimum between the
663%    main and the embedded scheme orders plus one
664ros_ELO = 2.0;
665%~~~> Y_stage_i ~ Y( T + H*Alpha_i )
666ros_Alpha(1) = 0.0;
667ros_Alpha(2) = 1.0;
668%~~~> Gamma_i = \sum_j  gamma_{i,j}
669ros_Gamma(1) = g;
670ros_Gamma(2) =-g;
671
672params = { ros_S, rosMethod, ros_A, ros_C, ros_M, ros_E, ...
673    ros_Alpha, ros_Gamma, ros_ELO, ros_NewF, ros_Name };
674
675return
676
677
678%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
679function [ params ] = Ros3()
680%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
681% --- AN L-STABLE METHOD, 3 stages, order 3, 2 function evaluations
682%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
683
684rosMethod = 2;
685%~~~> Name of the method
686ros_Name = 'ROS-3';
687%~~~> Number of stages
688ros_S = 3;
689
690%~~~> The coefficient matrices A and C are strictly lower triangular.
691%   The lower triangular (subdiagonal) elements are stored in row-wise order:
692%   A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
693%   The general mapping formula is:
694%       A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
695%       C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
696
697ros_A(1)= 1.0;
698ros_A(2)= 1.0;
699ros_A(3)= 0.0;
700
701ros_C(1) = -0.10156171083877702091975600115545E+01;
702ros_C(2) =  0.40759956452537699824805835358067E+01;
703ros_C(3) =  0.92076794298330791242156818474003E+01;
704%~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
705%   or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
706ros_NewF(1) = true;
707ros_NewF(2) = true;
708ros_NewF(3) = false;
709%~~~> M_i = Coefficients for new step solution
710ros_M(1) =  0.1E+01;
711ros_M(2) =  0.61697947043828245592553615689730E+01;
712ros_M(3) = -0.42772256543218573326238373806514;
713% E_i = Coefficients for error estimator
714ros_E(1) =  0.5;
715ros_E(2) = -0.29079558716805469821718236208017E+01;
716ros_E(3) =  0.22354069897811569627360909276199;
717%~~~> ros_ELO = estimator of local order - the minimum between the
718%    main and the embedded scheme orders plus 1
719ros_ELO = 3.0;
720%~~~> Y_stage_i ~ Y( T + H*Alpha_i )
721ros_Alpha(1)= 0.0;
722ros_Alpha(2)= 0.43586652150845899941601945119356;
723ros_Alpha(3)= 0.43586652150845899941601945119356;
724%~~~> Gamma_i = \sum_j  gamma_{i,j}
725ros_Gamma(1)= 0.43586652150845899941601945119356;
726ros_Gamma(2)= 0.24291996454816804366592249683314;
727ros_Gamma(3)= 0.21851380027664058511513169485832E+01;
728
729params = { ros_S, rosMethod, ros_A, ros_C, ros_M, ros_E, ...
730    ros_Alpha, ros_Gamma, ros_ELO, ros_NewF, ros_Name };
731
732return
733
734%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
735
736
737%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
738function [ params ] = Ros4()
739%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
740%     L-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 4 STAGES
741%     L-STABLE EMBEDDED ROSENBROCK METHOD OF ORDER 3
742%
743%      E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
744%      EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
745%      SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
746%      SPRINGER-VERLAG (1990)
747%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
748
749rosMethod = 3;
750%~~~> Name of the method
751ros_Name = 'ROS-4';
752%~~~> Number of stages
753ros_S = 4;
754
755%~~~> The coefficient matrices A and C are strictly lower triangular.
756%   The lower triangular (subdiagonal) elements are stored in row-wise order:
757%   A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
758%   The general mapping formula is:
759%       A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
760%       C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
761
762ros_A(1) = 0.2000000000000000E+01;
763ros_A(2) = 0.1867943637803922E+01;
764ros_A(3) = 0.2344449711399156;
765ros_A(4) = ros_A(2);
766ros_A(5) = ros_A(3);
767ros_A(6) = 0.0;
768
769ros_C(1) =-0.7137615036412310E+01;
770ros_C(2) = 0.2580708087951457E+01;
771ros_C(3) = 0.6515950076447975;
772ros_C(4) =-0.2137148994382534E+01;
773ros_C(5) =-0.3214669691237626;
774ros_C(6) =-0.6949742501781779;
775%~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
776%   or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
777ros_NewF(1)  = true;
778ros_NewF(2)  = true;
779ros_NewF(3)  = true;
780ros_NewF(4)  = false;
781%~~~> M_i = Coefficients for new step solution
782ros_M(1) = 0.2255570073418735E+01;
783ros_M(2) = 0.2870493262186792;
784ros_M(3) = 0.4353179431840180;
785ros_M(4) = 0.1093502252409163E+01;
786%~~~> E_i  = Coefficients for error estimator
787ros_E(1) =-0.2815431932141155;
788ros_E(2) =-0.7276199124938920E-01;
789ros_E(3) =-0.1082196201495311;
790ros_E(4) =-0.1093502252409163E+01;
791%~~~> ros_ELO  = estimator of local order - the minimum between the
792%    main and the embedded scheme orders plus 1
793ros_ELO  = 4.0;
794%~~~> Y_stage_i ~ Y( T + H*Alpha_i )
795ros_Alpha(1) = 0.0;
796ros_Alpha(2) = 0.1145640000000000E+01;
797ros_Alpha(3) = 0.6552168638155900;
798ros_Alpha(4) = ros_Alpha(3);
799%~~~> Gamma_i = \sum_j  gamma_{i,j}
800ros_Gamma(1) = 0.5728200000000000;
801ros_Gamma(2) =-0.1769193891319233E+01;
802ros_Gamma(3) = 0.7592633437920482;
803ros_Gamma(4) =-0.1049021087100450;
804
805params = { ros_S, rosMethod, ros_A, ros_C, ros_M, ros_E, ...
806    ros_Alpha, ros_Gamma, ros_ELO, ros_NewF, ros_Name };
807
808return
809
810%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
811function [ params ] = Rodas3()
812%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
813% --- A STIFFLY-STABLE METHOD, 4 stages, order 3
814%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
815
816rosMethod = 4;
817%~~~> Name of the method
818ros_Name = 'RODAS-3';
819%~~~> Number of stages
820ros_S = 4;
821
822%~~~> The coefficient matrices A and C are strictly lower triangular.
823%   The lower triangular (subdiagonal) elements are stored in row-wise order:
824%   A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
825%   The general mapping formula is:
826%       A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
827%       C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
828
829ros_A(1) = 0.0;
830ros_A(2) = 2.0;
831ros_A(3) = 0.0;
832ros_A(4) = 2.0;
833ros_A(5) = 0.0;
834ros_A(6) = 1.0;
835
836ros_C(1) = 4.0;
837ros_C(2) = 1.0;
838ros_C(3) =-1.0;
839ros_C(4) = 1.0;
840ros_C(5) =-1.0;
841ros_C(6) =-(8.0/3.0);
842
843%~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
844%   or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
845ros_NewF(1)  = true;
846ros_NewF(2)  = false;
847ros_NewF(3)  = true;
848ros_NewF(4)  = true;
849%~~~> M_i = Coefficients for new step solution
850ros_M(1) = 2.0;
851ros_M(2) = 0.0;
852ros_M(3) = 1.0;
853ros_M(4) = 1.0;
854%~~~> E_i  = Coefficients for error estimator
855ros_E(1) = 0.0;
856ros_E(2) = 0.0;
857ros_E(3) = 0.0;
858ros_E(4) = 1.0;
859%~~~> ros_ELO  = estimator of local order - the minimum between the
860%    main and the embedded scheme orders plus 1
861ros_ELO  = 3.0;
862%~~~> Y_stage_i ~ Y( T + H*Alpha_i )
863ros_Alpha(1) = 0.0;
864ros_Alpha(2) = 0.0;
865ros_Alpha(3) = 1.0;
866ros_Alpha(4) = 1.0;
867%~~~> Gamma_i = \sum_j  gamma_{i,j}
868ros_Gamma(1) = 0.5;
869ros_Gamma(2) = 1.5;
870ros_Gamma(3) = 0.0;
871ros_Gamma(4) = 0.0;
872
873params = { ros_S, rosMethod, ros_A, ros_C, ros_M, ros_E, ...
874    ros_Alpha, ros_Gamma, ros_ELO, ros_NewF, ros_Name };
875
876return
877
878%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
879function [ params ] = Rodas4()
880%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
881%     STIFFLY-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 6 STAGES
882%
883%      E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
884%      EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
885%      SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
886%      SPRINGER-VERLAG (1996)
887%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
888
889rosMethod = 5;
890%~~~> Name of the method
891ros_Name = 'RODAS-4';
892%~~~> Number of stages
893ros_S = 6;
894
895%~~~> Y_stage_i ~ Y( T + H*Alpha_i )
896ros_Alpha(1) = 0.000;
897ros_Alpha(2) = 0.386;
898ros_Alpha(3) = 0.210;
899ros_Alpha(4) = 0.630;
900ros_Alpha(5) = 1.000;
901ros_Alpha(6) = 1.000;
902
903%~~~> Gamma_i = \sum_j  gamma_{i,j}
904ros_Gamma(1) = 0.2500000000000000;
905ros_Gamma(2) =-0.1043000000000000;
906ros_Gamma(3) = 0.1035000000000000;
907ros_Gamma(4) =-0.3620000000000023E-01;
908ros_Gamma(5) = 0.0;
909ros_Gamma(6) = 0.0;
910
911%~~~> The coefficient matrices A and C are strictly lower triangular.
912%   The lower triangular (subdiagonal) elements are stored in row-wise order:
913%   A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
914%   The general mapping formula is:  A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
915%                  C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
916
917ros_A(1) = 0.1544000000000000E+01;
918ros_A(2) = 0.9466785280815826;
919ros_A(3) = 0.2557011698983284;
920ros_A(4) = 0.3314825187068521E+01;
921ros_A(5) = 0.2896124015972201E+01;
922ros_A(6) = 0.9986419139977817;
923ros_A(7) = 0.1221224509226641E+01;
924ros_A(8) = 0.6019134481288629E+01;
925ros_A(9) = 0.1253708332932087E+02;
926ros_A(10) =-0.6878860361058950;
927ros_A(11) = ros_A(7);
928ros_A(12) = ros_A(8);
929ros_A(13) = ros_A(9);
930ros_A(14) = ros_A(10);
931ros_A(15) = 1.0;
932
933ros_C(1) =-0.5668800000000000E+01;
934ros_C(2) =-0.2430093356833875E+01;
935ros_C(3) =-0.2063599157091915;
936ros_C(4) =-0.1073529058151375;
937ros_C(5) =-0.9594562251023355E+01;
938ros_C(6) =-0.2047028614809616E+02;
939ros_C(7) = 0.7496443313967647E+01;
940ros_C(8) =-0.1024680431464352E+02;
941ros_C(9) =-0.3399990352819905E+02;
942ros_C(10) = 0.1170890893206160E+02;
943ros_C(11) = 0.8083246795921522E+01;
944ros_C(12) =-0.7981132988064893E+01;
945ros_C(13) =-0.3152159432874371E+02;
946ros_C(14) = 0.1631930543123136E+02;
947ros_C(15) =-0.6058818238834054E+01;
948
949%~~~> M_i = Coefficients for new step solution
950ros_M(1) = ros_A(7);
951ros_M(2) = ros_A(8);
952ros_M(3) = ros_A(9);
953ros_M(4) = ros_A(10);
954ros_M(5) = 1.0;
955ros_M(6) = 1.0;
956
957%~~~> E_i  = Coefficients for error estimator
958ros_E(1) = 0.0;
959ros_E(2) = 0.0;
960ros_E(3) = 0.0;
961ros_E(4) = 0.0;
962ros_E(5) = 0.0;
963ros_E(6) = 1.0;
964
965%~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
966%   or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
967ros_NewF(1) = true;
968ros_NewF(2) = true;
969ros_NewF(3) = true;
970ros_NewF(4) = true;
971ros_NewF(5) = true;
972ros_NewF(6) = true;
973
974%~~~> ros_ELO  = estimator of local order - the minimum between the
975%        main and the embedded scheme orders plus 1
976ros_ELO = 4.0;
977
978params = { ros_S, rosMethod, ros_A, ros_C, ros_M, ros_E, ...
979    ros_Alpha, ros_Gamma, ros_ELO, ros_NewF, ros_Name };
980
981return
982
983% End of INTEGRATE function
984% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
985
986
Note: See TracBrowser for help on using the repository browser.