Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int/runge_kutta.m @ 4127

Last change on this file since 4127 was 2696, checked in by kanani, 7 years ago

Merge of branch palm4u into trunk
File size: 48.1 KB
Line 
1function [T, Y, RCNTRL, ICNTRL, RSTAT, ISTAT] = ...
2    RK_Int(Function, Tspan, Y0, Options, RCNTRL, ICNTRL)
3%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
4%  Implementation of Fully Implicit RK methods with the
5%  Coefficients:
6%               * Radau2A                                               
7%               * Lobatto3C                                             
8%               * Radau1A                                               
9%               * Gauss                                             
10%               * Lobatto3A                                             
11%                                                                     
12%    Solves the system y'=F(t,y) using a Fully Implicit RK method.
13%
14%    For details on Fully Implicit RK methods and their implementation consult:
15%      E. Hairer and G. Wanner
16%      "Solving ODEs II. Stiff and differential-algebraic problems".
17%      Springer series in computational mathematics, Springer-Verlag, 1996.
18%    The codes contained in the book inspired this implementation.
19%
20%    MATLAB implementation (C) Vishwas Rao (visrao@vt.edu).     
21%    Virginia Polytechnic Institute and State University             
22%    March, 2011   
23%
24%    Based on the Fortran90 implementation (C) Adrian Sandu, August 2004
25%    and revised by Philipp Miehe and Adrian Sandu, May 2006.
26%    Virginia Polytechnic Institute and State University
27%    Contact: sandu@cs.vt.edu
28%
29%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
30%  Input Arguments :
31%  The first four arguments are similar to the input arguments of
32%  MATLAB's ODE solvers
33%      Function  - A function handle for the ODE function
34%      Tspan     - The time space to integrate
35%      Y0        - Initial value
36%      Options   - ODE solver options created by odeset():           
37%                  AbsTol, InitialStep, Jacobian, MaxStep, and RelTol
38%                  are considered.  Other options are ignored.       
39%                  'Jacobian' must be set.                           
40%      RCNTRL    - real value input parameters (explained below)
41%      ICNTRL    - integer input parameters (explained below)
42%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
43%  Output Arguments:
44%  The first two arguments are similar to the output arguments of
45%  MATLAB's ODE solvers
46%      T         - A vector of final integration times.
47%      Y         - A matrix of function values.  Y(T(i),:) is the value of
48%                  the function at the ith output time.
49%      RCNTRL    - real value input parameters (explained below)
50%      ICNTRL    - integer input parameters (explained below)
51%      RSTAT     - real output parameters (explained below)
52%      ISTAT     - integer output parameters (explained below)
53%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
54%
55%    RCNTRL and ICNTRL on input and output:
56%
57%    Note: For input parameters equal to zero the default values of the
58%       corresponding variables are used.
59%
60%   
61%    ICNTRL(1) = 0: AbsTol, RelTol are N-dimensional vectors
62%              = 1: AbsTol, RelTol are scalars
63%
64%    ICNTRL(2)  -> selection of coefficients
65%        = 0 :    Radau2A
66%        = 1 :    Radau2A
67%        = 2 :    Lobatto3C
68%        = 3 :    Gauss
69%        = 4 :    Radau1A
70%        = 5 :    Lobatto3A
71%       
72%       
73%
74%    ICNTRL(4)  -> maximum number of integration steps
75%        For ICNTRL(4)=0) the default value of 100000 is used
76%
77%    RCNTRL(1)  -> Hmin, lower bound for the integration step size
78%          It is strongly recommended to keep Hmin = ZERO
79%    RCNTRL(2)  -> Hmax, upper bound for the integration step size
80%    RCNTRL(3)  -> Hstart, starting value for the integration step size
81%
82%    RCNTRL(4)  -> FacMin, lower bound on step decrease factor (default=0.2)
83%    RCNTRL(5)  -> FacMax, upper bound on step increase factor (default=6)
84%    RCNTRL(6)  -> FacRej, step decrease factor after multiple rejections
85%                          (default=0.1)
86%    RCNTRL(7)  -> FacSafe, by which the new step is slightly smaller
87%         than the predicted value  (default=0.9)
88%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
89%
90%    RSTAT and ISTAT on output:
91%
92%    ISTAT(1)  -> No. of function calls
93%    ISTAT(2)  -> No. of jacobian calls
94%    ISTAT(3)  -> No. of steps
95%    ISTAT(4)  -> No. of accepted steps
96%    ISTAT(5)  -> No. of rejected steps (except at very beginning)
97%    ISTAT(6)  -> No. of LU decompositions
98%    ISTAT(7)  -> No. of forward/backward substitutions
99%    ISTAT(8)  -> No. of singular matrix decompositions
100%
101%    RSTAT(1)  -> Texit, the time corresponding to the
102%                     computed Y upon return
103%    RSTAT(2)  -> Hexit, last accepted step before exit
104%    RSTAT(3)  -> Hnew, last predicted step (not yet taken)
105%                   For multiple restarts, use Hnew as Hstart
106%                     in the subsequent run
107%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
108global Nfun Njac Nstp Nacc Nrej Ndec Nsol Nsng Ntexit Nhexit Nhnew
109% Parse ODE options
110AbsTol = odeget(Options, 'AbsTol');
111if isempty(AbsTol)
112    AbsTol = 1.0e-3;
113end
114Hstart = odeget(Options, 'InitialStep');
115if isempty(Hstart)
116    Hstart = 0;
117end
118Jacobian = odeget(Options, 'Jacobian');
119if isempty(Jacobian)
120    error('A Jacobian function is required.');
121end
122Hmax = odeget(Options, 'MaxStep');
123if isempty(Hmax)
124    Hmax = 0;
125end
126RelTol = odeget(Options, 'RelTol');
127if isempty(RelTol)
128    RelTol = 1.0e-4;
129end
130% Initialize statistics
131Nfun=1; Njac=2; Nstp=3; Nacc=4; Nrej=5; Ndec=6; Nsol=7; Nsng=8;
132Ntexit=1; Nhexit=2; Nhnew=3;
133RSTATUS = zeros(20,1);
134ISTATUS = zeros(20,1);
135
136% Get problem size
137steps = length(Tspan);
138N = max(size(Y0));
139
140% Initialize tolerances
141ATOL = ones(N,1)*AbsTol;
142RTOL = ones(N,1)*RelTol;
143
144% Initialize step
145RCNTRL(2) = max(0, Hmax);
146RCNTRL(3) = max(0, Hstart);
147
148% Integrate
149Y = zeros(N,steps);
150T = zeros(steps,1);
151Y(:,1) = Y0;
152T(1) = Tspan(1);
153t=cputime;
154for i=2:steps
155    [T(i), Y(:,i), IERR,ISTATUS,RSTATUS] = ...
156        RK_Integrate(N, Y(:,i-1), T(i-1), Tspan(i), ...
157            Function, Jacobian, ATOL, RTOL, RCNTRL, ICNTRL,ISTATUS,RSTATUS);
158
159    if IERR < 0
160        error(['IRK exited with IERR=',num2str(IERR)]);
161    end
162
163end
164cputime-t
165RSTAT=RSTATUS;
166ISTAT=ISTATUS;
167%~~~>  Statistics on the work performed by the RK method
168function [T,Y,Ierr,ISTATUS,RSTATUS] = RK_Integrate(N,Y0,TIN, TOUT,OdeFunction1,...
169    OdeJacobian1,ATOL,RTOL,RCNTRL,ICNTRL,ISTATUS,RSTATUS)
170global ZERO ONE Nfun Njac Nstp Nacc Nrej Ndec Nsol Nsng Ntexit Nhacc Nhnew
171global R2A R1A L3C GAU L3A RKmax
172ZERO =0;ONE=1;Nfun=1;Njac=2;Nstp=3;Nacc=4;Nrej=5;Ndec=6;Nsol=7;Nsng=8;
173Ntexit=1;Nhacc=2;Nhnew=3;RKmax=3;R2A=1;R1A=2;L3C=3;GAU=4;L3A=5;
174Ierr =0;
175RSTATUS = zeros(20,1);
176ISTATUS = zeros(20,1);
177T1 = TIN;
178T2 = TOUT;
179Y(:,1)=Y0;
180[T,Y,Ierr,RSTATUS,ISTATUS] = RungeKutta(N, T1,T2,Y0,RTOL,ATOL,RCNTRL,ICNTRL,RSTATUS,ISTATUS,Ierr,OdeFunction1,OdeJacobian1);
181
182
183
184
185
186if(Ierr < 0)
187    disp('Runge-kutta Unsuccessful exit at T=',num2str(TIN));
188end
189return
190
191function [T,Y,Ierr,RSTATUS,ISTATUS]=RungeKutta(N,T1,T2,Y,RelTol,AbsTol,RCNTRL,...
192    ICNTRL,ISTATUS,RSTATUS,Ierr,OdeFunction1,OdeJacobian1)
193global ZERO  Roundoff SdirkError
194T=T1;
195
196if (ICNTRL(1)==0)
197            ITOL=1;
198else
199            ITOL=0;
200end
201
202if(ICNTRL(9)==0)
203    SdirkError = 0;
204else
205    SdirkError = 1;
206end
207
208switch(ICNTRL(2))
209        case{0,1}
210                Radau2A_Coefficients();
211        case(2)
212                Lobatto3C_Coefficients();
213        case(3)
214                Gauss_Coefficients();
215        case(4)
216                Radau1A_Coefficients();
217        case(5)
218                Lobatto3A_Coefficients();
219        otherwise
220                disp(['ICNTRL(2)=',num2str(ICNTRL(2))]);
221                RK_ErrorMsg(-13,Tin,0,IERR);
222                return;
223end
224
225if (ICNTRL(3)==0)
226                Max_no_steps=200000;
227else
228                Max_no_steps=ICNTRL(3);
229        if(Max_no_steps<0)
230            disp(['User-selected max no. of steps is -ve and is ',num2str(ICNTRL(3))]);
231        end
232end
233
234if (ICNTRL(4)==0)
235        NewtonMaxit=8;
236else
237        NewtonMaxit=ICNTRL(4);
238        if NewtonMaxit<=0
239                disp(['User-selected max no. of newton iterations is -ve and is ',num2str(ICNTRL(5))]);
240                RK_ErrorMsg(-2, T, ZERO, IERR);
241        end
242end
243
244if (ICNTRL(5)==0)
245    StartNewton=1;
246else
247    StartNewton =0;
248end
249
250if(ICNTRL(10)==0)
251    Gustafsson = 1;
252else
253    Gustafsson = 0;
254end
255Roundoff = eps;
256if(RCNTRL(1)==ZERO)
257    Hmin = ZERO;
258else
259    Hmin = min(abs(RCNTRL(1)),abs(T2 -T1));
260end
261if(RCNTRL(2)==ZERO)
262    Hmax=abs(T2-T1);
263else
264        Hmax= min(abs(RCNTRL(2)),abs(T2-T1));
265end
266if RCNTRL(3)==0
267        Hstart=0;
268else
269        Hstart = min(abs(RCNTRL(3)),abs(T2-T1));
270end
271
272%FacMin-- Lower bound on step decrease factor
273if (RCNTRL(4) == 0)
274        FacMin = 0.2;
275else
276        FacMin = RCNTRL(4);
277end
278%FacMax--Upper bound on step increase factor
279if RCNTRL(5)==0
280        FacMax=8.0;
281else
282        FacMax=RCNTRL(5);
283end
284%FacRej--step decrease factor after 2 consecutive rejections
285if RCNTRL(6)==0
286        FacRej=0.1;
287else
288        FacRej=RCNTRL(6);
289end
290
291%Facsafe:by which the new step is slightly smaller than the
292        %predicted value
293if RCNTRL(7)==0
294        FacSafe=0.9;
295else
296        FacSafe=RCNTRL(7);
297end
298
299if ( (FacMax < 1) ||( FacMin > 1) || (FacSafe <= 1.0e-03) || (FacSafe >= 1) )
300        disp(['\n RCNTRL(5)=',num2str(RCNTRL(4)) ' RCNTRL[5]=',num2str(RCNTRL(5)) ' RCNTRL[6]=', num2str(RCNTRL(6)) ' RCNTRL[7]=',num2str(RCNTRL(7))]);
301        RK_ErrorMsg(-4, T, ZERO, Ierr);
302end
303
304%ThetaMin: Decides whether the jacobian should be recomputed
305if RCNTRL(8)==0
306        ThetaMin=1.0e-03;
307else
308        ThetaMin=RCNTRL(8);
309end
310if (ThetaMin <= 0.0 || ThetaMin >= 1.0)
311        disp(['RCNTRL[8]=', num2str(RCNTRL(8))]);
312        RK_ErrorMsg(-5, Tin, ZERO, Ierr);
313end
314
315if RCNTRL(9)==0
316        NewtonTol = 3.0e-02;
317else
318        NewtonTol = RCNTRL(9);
319        if NewtonTol<=Roundoff
320                disp(['RCNTRL(9)=',num2str(RCNTRL(9))]);
321                RK_ErrorMsg(-6, Tin, ZERO,Ierr);
322        end
323end
324
325%Qmin and Qmax: If Qmin < Hnew/Hold < Qmax then step size is constant
326if RCNTRL(10)==0
327        Qmin=1;
328else
329        Qmin = RCNTRL(10);
330end
331if RCNTRL(11)==0
332        Qmax=1.2;
333else
334        Qmax=RCNTRL(11);
335end
336if (Qmin > 1 || Qmax <1)
337        disp(['Qmin=', num2str(RCNTRL(10)) ' Qmax=',num2str(RCNTRl(11))]);
338        RK_ErrorMsg(-7, T, ZERO, Ierr);
339end
340
341%check if tolerances are reasonable
342if ITOL==0
343        if AbsTol(1)<=0 ||  RelTol(1) <= 10.0*Roundoff
344                disp(['AbsTol=', num2str(AbsTol(1)) ' RelTol=',num2str(RelTol(1))]);
345                RK_ErrorMsg(-8, T, ZERO, Ierr);
346        end
347else
348        for i=1:N
349        if (( AbsTol(i) <= 0) || (RelTol(i) <= 10.0*Roundoff) )
350            disp(['AbsTol(',num2str(i) ')= ',num2str(AbsTol(i))]);
351                        disp(['RelTol(',num2str(i) ')=  ',RelTol(i)]);
352                        RK_ErrorMsg(-8, T, ZERO, Ierr);
353        end
354        end
355end
356
357if(Ierr < 0)
358    return;
359end
360%Call the core Integrator
361[T,Y,Ierr,ISTATUS,RSTATUS] = RK_Integrator(N,T1,T2,Y,Hstart,Hmin,...
362                            Hmax,Roundoff,AbsTol,RelTol,...
363                            ITOL,Max_no_steps,StartNewton,NewtonTol,ThetaMin,...
364                            FacSafe,FacMin,FacMax,FacRej,Qmin,Qmax,NewtonMaxit,ISTATUS,RSTATUS,Gustafsson,Ierr,OdeFunction1,OdeJacobian1);
365return;
366
367
368
369function [T,Y,Ierr,ISTATUS,RSTATUS] = RK_Integrator(N,T1,T2,Y,Hstart,Hmin,Hmax,...
370                                      Roundoff,AbsTol,RelTol,ITOL,Max_no_steps,...
371                                      StartNewton,NewtonTol,ThetaMin,FacSafe,FacMin,...
372                                      FacMax,FacRej,Qmin,Qmax,NewtonMaxit,ISTATUS,...
373                                      RSTATUS,Gustafsson,Ierr,OdeFunction1,OdeJacobian1)
374global SdirkError ZERO Nfun Njac Nsng  Nstp ONE rkMethod L3A rkBgam
375global rkTheta rkGamma rkD rkELO Nacc Ntexit Nhacc Nhnew Hacc
376T = T1;
377CONT=zeros(N,3);
378Tdirection = sign(T2-T);
379H = min (max(abs(Hmin),abs(Hstart)),Hmax);
380if(abs(H)<=10*Roundoff)
381    H=1.0e-6;
382end
383H =sign(Tdirection)*abs(H);
384Hold = H;
385Reject = 0;
386FirstStep = 1;
387SkipJac = 0;
388SkipLU = 0;
389if((T+H*1.0001-T2)*Tdirection>=ZERO)
390    H = T2 -T;
391end
392Nconsecutive = 0;
393SCAL = RK_ErrorScale(N,ITOL,AbsTol,RelTol,Y);
394while((T2-T)*Tdirection-Roundoff>ZERO)
395    F0 = OdeFunction1(T,Y);
396    ISTATUS(Nfun)=ISTATUS(Nfun)+1;
397    if(SkipLU == 0)
398        if(SkipJac==0)
399            FJAC = OdeJacobian1(T,Y);
400            ISTATUS(Njac)=ISTATUS(Njac)+1;
401        end
402        [E1L,E1U,E2L,E2U,ISING,ISTATUS]=RK_Decomp(N,H,FJAC,ISTATUS);
403        if(ISING ~= 0)
404            ISTATUS(Nsng)=ISTATUS(Nsng)+1;
405            Nconsecutive=Nconsecutive+1;
406            if(Nconsecutive>=5)
407                RK_ErrorMsg(-12,T,H,Ierr);
408            end
409            H=H*0.5;
410            Reject = 1;
411            SkipJac = 1;
412            SkipLU = 0;
413            continue;
414        else
415            Nconsecutive = 0;
416        end
417    end
418    ISTATUS(Nstp)=ISTATUS(Nstp)+1;
419    if (ISTATUS(Nstp)>Max_no_steps)
420                disp(['Max number of time steps is =', num2str(Max_no_steps)]);
421                RK_ErrorMsg(-9,T,H,IERR);
422    end
423    if (0.1*abs(H)<=abs(T)*Roundoff)
424                RK_ErrorMsg(-10,T,H,Ierr);
425    end
426    %Loop for simplified newton iterations
427
428        %Starting values for newton iterations
429    if((FirstStep==1)||(StartNewton==0))
430                        Z1=zeros(N,1);
431                        Z2=zeros(N,1);
432                        Z3=zeros(N,1);
433        else
434                evaluate=2;
435                [CONT, Z1, Z2, Z3]=RK_Interpolate(evaluate, N, H, Hold,Z1,Z2,Z3,CONT);
436    end
437   
438    %/*~~~> Initializations for Newton iteration */
439   
440    NewtonDone = 0;
441    Fac = 0.5;
442    for NewtonIter=1:NewtonMaxit
443        %Preparing RHS
444        [R1,R2,R3]=RK_PrepareRHS(N,T,H,Y,F0,Z1,Z2,Z3,OdeFunction1);
445        %Solve the linear systems
446        [R1,R2,R3,ISTATUS] = RK_Solve(N,H,E1L,E1U,E2L,E2U,R1,R2,R3,ISTATUS);
447        %Find NewtonIncrement
448        NewtonIncrement=sqrt((RK_ErrorNorm(N,SCAL,R1)*RK_ErrorNorm(N,SCAL,R1)...
449                                                                        +RK_ErrorNorm(N,SCAL,R2)*RK_ErrorNorm(N,SCAL,R2)...
450                                                                        +RK_ErrorNorm(N,SCAL,R3)*RK_ErrorNorm(N,SCAL,R3))...
451                                                        /3.0);
452        if(NewtonIter == 1)
453            Theta = abs(ThetaMin);
454            NewtonRate = 2.0;
455        else
456            Theta = NewtonIncrement/NewtonIncrementOld;
457            if(Theta < 0.99)
458                NewtonRate = Theta/(ONE-Theta);
459            else
460                break
461            end
462            if(NewtonIter<NewtonMaxit)
463                NewtonPredictedErr = NewtonIncrement*Theta^(NewtonMaxit-NewtonIter)/(ONE - Theta);
464                if(NewtonPredictedErr>=NewtonTol)
465                    Qnewton = min(10.0,NewtonPredictedErr/NewtonTol);
466                    Fac = 0.8*Qnewton^(-ONE/(1+NewtonMaxit-NewtonIter));
467                    break
468                end
469            end
470        end
471        NewtonIncrementOld=max(NewtonIncrement,Roundoff);
472        Z1=Z1-R1;
473        Z2=Z2-R2;
474        Z3=Z3-R3;
475        NewtonDone = (NewtonRate*NewtonIncrement<=NewtonTol);
476        if(NewtonDone==1)
477            break;
478        end
479        if(NewtonIter==NewtonMaxit)
480            disp('Slow or no convergence in newton iterations');
481                        disp('Max no of newton iterations reached');
482        end
483    end
484    if(NewtonDone == 0)
485        H=Fac*H;
486        Reject=1;
487        SkipJac=1;
488                SkipLU=0;
489        continue;
490    end
491    %/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
492    %/*~~~> SDIRK Stage */
493    %/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
494    if(SdirkError==1)
495        Z4=Z3;
496        G=zeros(N,1);
497        if(rkMethod ~= L3A)
498            G=G+rkBgam(1)*H*F0;
499        end
500        G=G+Z1*rkTheta(1);
501        G=G+Z2*rkTheta(2);
502        G=G+Z3*rkTheta(3);
503        NewtonDone = 0;
504        Fac = 0.5;
505        for NewtonIter=1:NewtonMaxit
506            TMP=Y+Z4;
507            R4=OdeFunction1(T+H,TMP);
508            ISTATUS(Nfun)=ISTATUS(Nfun)+1;
509            R4=R4-(rkGamma/H)*Z4;
510            R4=R4+(rkGamma/H)*G;
511            [R4,ISTATUS] = RK_KppSolve(E1L,E1U,R4,ISTATUS);
512            NewtonIncrement = RK_ErrorNorm(N,SCAL,R4);
513            if(NewtonIter==1)
514                ThetaSD=abs(ThetaMin);
515                NewtonRate = 2.0;
516            else
517                ThetaSD = NewtonIncrement/NewtonIncrementOld;
518                if(ThetaSD<0.99)
519                    NewtonRate = ThetaSD/(ONE - ThetaSD);
520                    NewtonPredictedErr=NewtonIncrement*ThetaSD^((NewtonMaxit-NewtonIter)/(1-ThetaSD));
521                    if(NewtonPredictedErr>=NewtonTol)
522                                                Qnewton=min(10.0,NewtonPredictedErr/NewtonTol);
523                                                Fac=0.8*Qnewton^(-1/(1+NewtonMaxit-NewtonIter));
524                                                break;
525                    end
526                else
527                    break;
528                end
529            end
530            NewtonIncrementOld = NewtonIncrement;
531            Z4 = Z4+R4;
532            %Check error in NewtonIterations
533            NewtonDone=(NewtonRate*NewtonIncrement<=NewtonTol);
534            if(NewtonDone==1)
535                break;
536            end
537        end
538        if(NewtonDone == 0)
539            H=Fac*H;
540            Reject =1;
541            SkipJac=1;
542            SkipLU = 0;
543            continue;
544        end
545    end
546    %/*~~~> End of Simplified SDIRK Newton iterations */
547
548    %/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
549    %/*~~~> Error estimation */
550    %/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
551    if(SdirkError==1)
552        R4=zeros(N,1);
553        if(rkMethod==L3A)
554            R4=H*rkF(1)*F0;
555            if(rkF(2)~=ZERO)
556                R4=R4+Z1*rkF(1);
557            end
558           
559            if(rkF(3)~=ZERO)
560                R4=R4+Z2*rkF(2);
561            end
562           
563            if(rkF(4)~=ZERO)
564                R4=R4+Z3*rkF(3);
565            end
566           
567            TMP = Y+Z4;
568            R1 = OdeFunction1(T+H,TMP);
569            R4=R4+H*rkBgam(5)*R1;
570        else
571            if(rkD(1)~=0)
572                R4=R4+rkD(1)*Z1;
573            end
574            if(rkD(2)~=ZERO)
575                R4=R4+rkD(2)*Z2;
576            end
577            if(rkD(3)~=ZERO)
578                R4=R4+rkD(3)*Z3;
579            end
580            R4=R4-Z4;
581        end
582                Err=RK_ErrorNorm(N,SCAL,R4);
583    else
584        [Err,ISTATUS]=RK_ErrorEstimate(N,H,T,Y,F0,E1L,E1U,Z1,Z2,Z3,SCAL,...
585                                                        FirstStep,Reject,ISTATUS);
586    end
587    Fac=Err^(-1/rkELO)*min(FacSafe,(1+2*NewtonMaxit)/(NewtonIter+2*NewtonMaxit));
588        Fac=min(FacMax,max(FacMin,Fac));
589        Hnew=Fac*H;
590    %/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
591    %/*~~~> Accept/reject step */
592    %/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
593   
594    %Accept
595    if(Err < ONE)
596        FirstStep =0;
597        ISTATUS(Nacc)=ISTATUS(Nacc)+1;
598        if(Gustafsson==ONE)
599            if(ISTATUS(Nacc)>1)
600                FacGus=FacSafe*(H/Hacc)*(Err*Err/ErrOld)^(-0.25);
601                                FacGus=min(FacMax,max(FacMin,FacGus));
602                                Fac=min(Fac,FacGus);
603                                Hnew=Fac*H;
604            end
605            Hacc=H;
606                        ErrOld=max(1.0e-02,Err);
607        end
608        Hold=H;
609                T=T+H;
610                if(rkD(1)~=0)
611                        Y=Y+Z1*rkD(1);
612                end
613                if(rkD(2)~=0)
614                        Y=Y+Z2*rkD(2);
615                end;
616                if(rkD(3)~=0)
617                        Y=Y+Z3*rkD(3);
618                end;
619        %Construct the solution quadratic interpolant Q(c_i)=Z_i,i=1:3
620                if(StartNewton==1)
621                        evaluate=1;
622                        [CONT, Z1, Z2, Z3]=RK_Interpolate(evaluate, N, H, Hold,Z1,Z2,Z3,CONT);
623                end
624                SCAL=RK_ErrorScale(N,ITOL,AbsTol,RelTol,Y);
625                RSTATUS(Ntexit)=T;
626                RSTATUS(Nhnew)=Hnew;
627                RSTATUS(Nhacc)=H;
628                Hnew=Tdirection*min(max(abs(Hnew),Hmin),Hmax);
629        if(Reject==1)
630                        Hnew=Tdirection*min(abs(Hnew),abs(H));
631        end;
632                Reject=0;
633        if((T + Hnew/Qmin - T2)*Tdirection>=0)
634                        H=T2-T;
635                else
636                        Hratio=Hnew/H;
637                        SkipLU=((Theta<=ThetaMin) && (Hratio>=Qmin) && (Hratio<=Qmax));
638                        if(SkipLU==0)
639                                H=Hnew;
640                        end;
641           
642        end;
643         SkipJac=0;
644               
645        else
646                        if((FirstStep==1) || (Reject==1))
647                                H=FacRej*H;
648                        else
649                                H=Hnew;
650                        end;
651                        Reject=1;
652                        SkipJac=1;
653                        SkipLU=0;
654                        if(ISTATUS(Nacc)>=1)
655                                ISTATUS(Nrej)=ISTATUS(Nrej)+1;
656                        end;
657                end;
658        end;
659        Ierr=1;
660
661return;
662
663function [Err,ISTATUS]=RK_ErrorEstimate(N,H,T,Y,F0,E1L,E1U,Z1,Z2,Z3,SCAL,...
664                                                        FirstStep,Reject,ISTATUS)
665global rkE rkMethod R1A GAU L3A
666HrkE1=rkE(2)/H;
667HrkE2=rkE(3)/H;
668HrkE3=rkE(4)/H;
669F2=HrkE1*Z1+HrkE2*Z2+HrkE3*Z3;
670TMP = rkE(1)*F0+F2;
671[TMP,ISTATUS]=RK_KppSolve(E1L,E1U,TMP,ISTATUS);
672if((rkMethod == R1A) || (rkMethod == GAU) || (rkMethod == L3A))
673        [TMP,ISTATUS]=RK_KppSolve(E1L,E1U,TMP,ISTATUS);
674end
675
676if (rkMethod == GAU)
677        [TMP,ISTATUS]=RK_KppSolve(E1L,E1U,TMP,ISTATUS);
678end
679Err=RK_ErrorNorm(N,SCAL,TMP);
680if(Err < 1)
681    return;
682end
683if((FirstStep==1 )||(Reject==1))
684    TMP=TMP+Y;
685    F1=OdeFunction1(T,TMP);
686    ISTATUS(Nfun)=ISTATUS(Nfun)+1;
687    TMP=F1+F2;
688    [TMP,ISTATUS]=RK_KppSolve(E1L,E1U,TMP,ISTATUS);
689    Err=RK_ErrorNorm(N,SCAL,TMP);
690end
691return
692
693function RK_ErrorMsg(Code, T, H, IERR)
694       
695        Code=IERR;
696        disp('Forced to exit from RungeKutta due to the following error:');
697        switch(Code)
698                case(-1)
699                        disp('--> Improper value for maximal no of steps');
700                case (-2)
701                disp('--> Improper value for maximal no of Newton iterations');
702                case (-3)
703                disp('--> Hmin/Hmax/Hstart must be positive');
704                case (-4)
705                disp('--> Improper values for FacMin/FacMax/FacSafe/FacRej');
706                case (-5)
707                disp('--> Improper value for ThetaMin');
708                case (-6)
709                disp('--> Newton stopping tolerance too small');
710                case (-7)
711                disp('--> Improper values for Qmin, Qmax');
712                case (-8)
713                disp('--> Tolerances are too small');
714                case (-9)
715                disp('--> No of steps exceeds maximum bound');
716                case (-10)
717                disp('--> Step size too small: (T + 10*H = T) or H < Roundoff');
718                case (-11)
719                disp('--> Matrix is repeatedly singular');
720                case (-12)
721                disp('--> Non-convergence of Newton iterations');
722                case (-13)
723                disp('--> Requested RK method not implemented');
724                otherwise
725                disp(['Unknown Error code:', num2str(Code)]);
726        end;
727
728        disp(['T=',num2str(T) 'H=',num2str(H)]);
729return;
730
731function SCAL=RK_ErrorScale(N,ITOL,AbsTol,RelTol,Y)
732        if(ITOL==0)
733        SCAL=1./(AbsTol(1)+RelTol(1)*abs(Y));
734    else
735        SCAL=1./(AbsTol+RelTol.*abs(Y));
736        end;
737return;
738
739function ErrorNorm = RK_ErrorNorm(N,SCAL,DY)
740
741    temp=DY.*DY.*SCAL.*SCAL;
742    ErrorNorm = sum(temp);
743ErrorNorm=max(sqrt(ErrorNorm/N),1.0e-10);
744
745function [CONT, Z1, Z2, Z3]=RK_Interpolate(action, N, H, Hold,Z1,Z2,Z3,CONT)
746        global rkC
747   
748
749        if(action==1) %Make
750                den=(rkC(3)-rkC(2))*(rkC(2)-rkC(1))*(rkC(1)-rkC(3));
751                for i=1:N
752                        CONT(i,1)=(-rkC(3)*rkC(3)*rkC(2)*Z1(i)...
753                                          +Z3(i)*rkC(2)*rkC(1)*rkC(1)...
754                                          +rkC(2)*rkC(2)*rkC(3)*Z1(i)...
755                                          -rkC(2)*rkC(2)*rkC(1)*Z3(i)...
756                                          +rkC(3)*rkC(3)*rkC(1)*Z2(i)...
757                                  -Z2(i)*rkC(3)*rkC(1)*rkC(1))/den-Z3(i);
758                        CONT(i,2) = -( rkC(1)*rkC(1)*(Z3(i)-Z2(i))...
759                        + rkC(2)*rkC(2)*(Z1(i)-Z3(i))...
760                        + rkC(3)*rkC(3)*(Z2(i)-Z1(i)) )/den;
761
762                        CONT(i,3) = ( rkC(1)*(Z3(i)-Z2(i))...
763                        + rkC(2)*(Z1(i)-Z3(i))...
764                        + rkC(3)*(Z2(i)-Z1(i)) )/den;
765                end;
766    end
767   
768    if (action==2) %Eval
769                r=H/Hold;
770                x1=1+rkC(1)*r;
771                x2 = 1 + rkC(2)*r;
772                x3 = 1 + rkC(3)*r;
773                for i=1:N
774               
775                        Z1(i) = CONT(i,1)+x1*(CONT(i,2)+x1*CONT(i,3)); 
776                        Z2(i) = CONT(i,1)+x2*(CONT(i,2)+x2*CONT(i,3)); 
777                        Z3(i) = CONT(i,1)+x3*(CONT(i,2)+x3*CONT(i,3));
778                end;
779    end;
780return;
781
782
783
784function [R1,R2,R3]=RK_PrepareRHS(N,T,H,Y,FO,Z1,Z2,Z3,OdeFunction1)
785        global rkMethod
786        global rkA rkC
787    global L3A
788        R1=Z1;
789        R2=Z2;
790        R3=Z3;
791        if(rkMethod==L3A)
792                R1=R1-H*rkA(1,1)*FO;
793                R2=R2-H*rkA(2,1)*FO;
794                R3=R3-H*rkA(3,1)*FO;
795        end
796        TMP=Y+Z1;
797        F=OdeFunction1(T+rkC(1)*H,TMP);
798        R1=R1-H*rkA(1,1)*F;
799        R2=R2-H*rkA(2,1)*F;
800        R3=R3-H*rkA(3,1)*F;
801
802        TMP=Y+Z2;
803        F=OdeFunction1(T+rkC(2)*H,TMP);
804        R1=R1-H*rkA(1,2)*F;
805        R2=R2-H*rkA(2,2)*F;
806        R3=R3-H*rkA(3,2)*F;
807
808        TMP=Y+Z3;
809        F=OdeFunction1(T+rkC(3)*H,TMP);
810        R1=R1-H*rkA(1,3)*F;
811        R2=R2-H*rkA(2,3)*F;
812        R3=R3-H*rkA(3,3)*F;
813return;
814
815function [E1L,E1U,E2L,E2U,ISING,ISTATUS]=RK_Decomp(N, H, FJAC,ISTATUS)
816        global Ndec
817        global rkGamma rkAlpha rkBeta
818    ISING =0;
819        Gamma = rkGamma / H;
820        Alpha = rkAlpha / H;
821        Beta  = rkBeta / H;
822        E1=Gamma*eye(N);
823    E1=E1-FJAC;
824
825
826       
827        [E1L,E1U]=lu(E1);
828    ISTATUS(Ndec)=ISTATUS(Ndec)+1;
829        if(det(E1L)==0 || det(E1U)==0)
830                ISING=1;
831        end;
832   
833        if(ISING ~= 0)
834               
835                return;
836        end;
837    E2R=complex(Alpha,Beta)*eye(N)-FJAC;
838       
839        [E2L,E2U]=lu(E2R);
840    ISTATUS(Ndec)=ISTATUS(Ndec)+1;
841        if(abs(det(E2L))==0 || abs(det(E2U))==0)
842                ISING=1;
843        end;
844        if(ISING ~= 0)
845                return;
846        end;
847   
848return
849
850function [R1,R2,R3,ISTATUS]=RK_Solve(N,H,E1L,E1U,E2L,E2U,R1,R2,R3,ISTATUS)
851        global Nsol 
852        global rkT rkTinvAinv
853        for i=1:N
854                x1=R1(i)/H;
855                x2=R2(i)/H;
856                x3=R3(i)/H;
857                R1(i) = rkTinvAinv(1,1)*x1 + rkTinvAinv(1,2)*x2 + rkTinvAinv(1,3)*x3;
858                R2(i) = rkTinvAinv(2,1)*x1 + rkTinvAinv(2,2)*x2 + rkTinvAinv(2,3)*x3;
859                R3(i) = rkTinvAinv(3,1)*x1 + rkTinvAinv(3,2)*x2 + rkTinvAinv(3,3)*x3;
860        end;
861        tmp = E1L\R1;
862        R1 =  E1U\tmp;
863        BCR=R2;
864        BCI=R3;
865        a1=complex(BCR,BCI);
866        tmp = E2L\a1;
867        BC = E2U\tmp;
868        R2=real(BC);
869        R3=imag(BC);
870        for i=1:N
871                x1 = R1(i);
872                x2 = R2(i);
873                x3 = R3(i);
874        R1(i) = rkT(1,1)*x1 + rkT(1,2)*x2 + rkT(1,3)*x3;
875                R2(i) = rkT(2,1)*x1 + rkT(2,2)*x2 + rkT(2,3)*x3;
876                R3(i) = rkT(3,1)*x1 + rkT(3,2)*x2 + rkT(3,3)*x3;
877        end
878        ISTATUS(Nsol)=ISTATUS(Nsol)+1;
879return
880
881function [R4,ISTATUS]=RK_KppSolve(E1L,E1U,R4,ISTATUS)
882global Nsol
883temp = E1L\R4;
884R4=E1U\temp;
885ISTATUS(Nsol)=ISTATUS(Nsol)+1;
886return
887
888
889
890
891function Radau2A_Coefficients ()
892
893        global rkMethod SdirkError
894        global rkT rkTinv rkTinvAinv rkAinvT rkA rkB rkC rkD rkE
895        global rkBgam  rkTheta  rkGamma rkAlpha rkBeta rkELO
896        global R2A
897
898       
899        if(SdirkError==1)
900                b0=0.2e-01;
901        else
902                b0=0.5e-01;
903        end;
904        rkMethod=R2A;
905        rkA(1,1) = 1.968154772236604258683861429918299e-01;
906        rkA(1,2) = -6.55354258501983881085227825696087e-02;
907        rkA(1,3) = 2.377097434822015242040823210718965e-02;
908        rkA(2,1) = 3.944243147390872769974116714584975e-01;
909        rkA(2,2) = 2.920734116652284630205027458970589e-01;
910        rkA(2,3) = -4.154875212599793019818600988496743e-02;
911        rkA(3,1) = 3.764030627004672750500754423692808e-01;
912        rkA(3,2) = 5.124858261884216138388134465196080e-01;
913        rkA(3,3) = 1.111111111111111111111111111111111e-01;
914
915        rkB(1) = 3.764030627004672750500754423692808e-01;
916        rkB(2) = 5.124858261884216138388134465196080e-01;
917        rkB(3) = 1.111111111111111111111111111111111e-01;
918
919        rkC(1) = 1.550510257216821901802715925294109e-01;
920        rkC(2) = 6.449489742783178098197284074705891e-01;
921        rkC(3) = 1;
922
923        % New solution: H* Sum B_j*f(Z_j) = Sum D_j*Z_j
924        rkD(1) =0;
925        rkD(2) =0;
926        rkD(3) =1;
927
928        % Classical error estimator: */
929        % H* Sum (B_j-Bhat_j)*f(Z_j) = H*E(0)*f(0) + Sum E_j*Z_j */
930        rkE(1) = 1*b0;
931        rkE(2) = -10.04880939982741556246032950764708*b0;
932        rkE(3) = 1.382142733160748895793662840980412*b0;
933        rkE(4) = -.3333333333333333333333333333333333*b0;
934
935    %   /* Sdirk error estimator */
936        rkBgam(1) = b0;
937        rkBgam(2) = .3764030627004672750500754423692807...
938                        -1.558078204724922382431975370686279*b0;
939        rkBgam(3) = 0.8914115380582557157653087040196118*b0...
940                        +0.5124858261884216138388134465196077;
941        rkBgam(4) = (-0.1637777184845662566367174924883037)...
942                        -0.3333333333333333333333333333333333*b0;
943        rkBgam(5) = 0.2748888295956773677478286035994148;
944
945        % H* Sum Bgam_j*f(Z_j) = H*Bgam(0)*f(0) + Sum Theta_j*Z_j */
946        rkTheta(1) = (-1.520677486405081647234271944611547)...
947                        -10.04880939982741556246032950764708*b0;
948        rkTheta(2) = (2.070455145596436382729929151810376)...
949                        +1.382142733160748895793662840980413*b0;
950        rkTheta(3) = -0.3333333333333333333333333333333333*b0...
951                        -.3744441479783868387391430179970741;
952
953        %/* Local order of error estimator */
954        if ( b0 == 0)
955                rkELO  = 6.0;
956        else
957                rkELO  = 4.0;
958        end;
959
960        %/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
961         % ~~~> Diagonalize the RK matrix:
962          %rkTinv * inv(rkA) * rkT =
963          %|  rkGamma      0           0     |
964          %|      0      rkAlpha   -rkBeta   |
965          %|      0      rkBeta     rkAlpha  |
966          %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
967
968        rkGamma = 3.637834252744495732208418513577775;
969        rkAlpha = 2.681082873627752133895790743211112;
970        rkBeta  = 3.050430199247410569426377624787569;
971
972        rkT(1,1) = 9.443876248897524148749007950641664e-02;
973        rkT(1,2) = -1.412552950209542084279903838077973e-01;
974        rkT(1,3) = -3.00291941051474244918611170890539e-02;
975        rkT(2,1) = 2.502131229653333113765090675125018e-01;
976        rkT(2,2) = 2.041293522937999319959908102983381e-01;
977        rkT(2,3) = 3.829421127572619377954382335998733e-01;
978        rkT(3,1) =  1;
979        rkT(3,2) =  1;
980        rkT(3,3) =  0;
981
982        rkTinv(1,1) = 4.178718591551904727346462658512057;
983        rkTinv(1,2) = 3.27682820761062387082533272429617e-01;
984        rkTinv(1,3) = 5.233764454994495480399309159089876e-01;
985        rkTinv(2,1) = -4.178718591551904727346462658512057;
986        rkTinv(2,2) = -3.27682820761062387082533272429617e-01;
987        rkTinv(2,3) = 4.766235545005504519600690840910124e-01;
988        rkTinv(3,1) = -5.02872634945786875951247343139544e-01;
989        rkTinv(3,2) = 2.571926949855605429186785353601676;
990        rkTinv(3,3) = -5.960392048282249249688219110993024e-01;
991
992        rkTinvAinv(1,1) = 1.520148562492775501049204957366528e+01;
993        rkTinvAinv(1,2) = 1.192055789400527921212348994770778;
994        rkTinvAinv(1,3) = 1.903956760517560343018332287285119;
995        rkTinvAinv(2,1) = -9.669512977505946748632625374449567;
996        rkTinvAinv(2,2) = -8.724028436822336183071773193986487;
997        rkTinvAinv(2,3) = 3.096043239482439656981667712714881;
998        rkTinvAinv(3,1) = -1.409513259499574544876303981551774e+01;
999        rkTinvAinv(3,2) = 5.895975725255405108079130152868952;
1000        rkTinvAinv(3,3) = -1.441236197545344702389881889085515e-01;
1001
1002        rkAinvT(1,1) = .3435525649691961614912493915818282;
1003        rkAinvT(1,2) = -.4703191128473198422370558694426832;
1004        rkAinvT(1,3) = .3503786597113668965366406634269080;
1005        rkAinvT(2,1) = .9102338692094599309122768354288852;
1006        rkAinvT(2,2) = 1.715425895757991796035292755937326;
1007        rkAinvT(2,3) = .4040171993145015239277111187301784;
1008        rkAinvT(3,1) = 3.637834252744495732208418513577775;
1009        rkAinvT(3,2) = 2.681082873627752133895790743211112;
1010        rkAinvT(3,3) = -3.050430199247410569426377624787569;
1011return;
1012
1013function Lobatto3C_Coefficients()
1014       
1015        global rkMethod SdirkError
1016        global rkT rkTinv rkTinvAinv rkAinvT rkA rkB rkC rkD rkE
1017        global rkBgam rkBhat rkTheta rkGamma rkAlpha rkBeta rkELO
1018        global L3C
1019
1020       
1021
1022        rkMethod=L3C;
1023        if(SdirkError==1)
1024                b0=0.2;
1025        else
1026                b0=0.5;
1027        end;
1028        rkA(1,1) = .1666666666666666666666666666666667;
1029        rkA(1,2) = -.3333333333333333333333333333333333;
1030        rkA(1,3) = .1666666666666666666666666666666667;
1031        rkA(2,1) = .1666666666666666666666666666666667;
1032        rkA(2,2) = .4166666666666666666666666666666667;
1033        rkA(2,3) = -.8333333333333333333333333333333333e-01;
1034        rkA(3,1) = .1666666666666666666666666666666667;
1035        rkA(3,2) = .6666666666666666666666666666666667;
1036        rkA(3,3) = .1666666666666666666666666666666667;
1037
1038        rkB(1) = .1666666666666666666666666666666667;
1039        rkB(2) = .6666666666666666666666666666666667;
1040        rkB(3) = .1666666666666666666666666666666667;
1041
1042        rkC(1) = 0;
1043        rkC(2) = 0.5;
1044        rkC(3) = 1;
1045
1046        %/* Classical error estimator, embedded solution: */
1047        rkBhat(1) = b0;
1048        rkBhat(2) = .16666666666666666666666666666666667-b0;
1049        rkBhat(3) = .66666666666666666666666666666666667;
1050        rkBhat(4) = .16666666666666666666666666666666667;
1051
1052%       /* New solution: h Sum_j b_j f(Z_j) = sum d_j Z_j */
1053        rkD(1) = 0;
1054        rkD(2) = 0;
1055        rkD(3) = 1;
1056
1057%       /* Classical error estimator: */
1058%       /* H* Sum (B_j-Bhat_j)*f(Z_j) = H*E(0)*f(0) + Sum E_j*Z_j */
1059        rkE(1) = .3808338772072650364017425226487022*b0;
1060        rkE(2) = (-1.142501631621795109205227567946107)*b0;
1061        rkE(3) = (-1.523335508829060145606970090594809)*b0;
1062        rkE(4) = .3808338772072650364017425226487022*b0;
1063
1064%       /* Sdirk error estimator */
1065        rkBgam(1) = b0;
1066        rkBgam(2) = .1666666666666666666666666666666667-1.0*b0;
1067        rkBgam(3) = .6666666666666666666666666666666667;
1068        rkBgam(4) = (-.2141672105405983697350758559820354);
1069        rkBgam(5) = .3808338772072650364017425226487021;
1070
1071%       /* H* Sum Bgam_j*f(Z_j) = H*Bgam(0)*f(0) + Sum Theta_j*Z_j */
1072        rkTheta(1) = -3.0*b0-.3808338772072650364017425226487021;
1073        rkTheta(2) = -4.0*b0+1.523335508829060145606970090594808;
1074        rkTheta(3) = (-.142501631621795109205227567946106)+b0;
1075
1076%       /* Local order of error estimator */
1077        if (b0 == 0)
1078                        rkELO  = 5.0;
1079        else
1080                        rkELO  = 4.0;
1081
1082        %/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1083        %  ~~~> Diagonalize the RK matrix:
1084        %  rkTinv * inv(rkA) * rkT =
1085        %  |  rkGamma      0           0     |
1086        %  |      0      rkAlpha   -rkBeta   |
1087        %  |      0      rkBeta     rkAlpha  |
1088        %  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1089
1090        rkGamma = 2.625816818958466716011888933765284;
1091        rkAlpha = 1.687091590520766641994055533117359;
1092        rkBeta  = 2.508731754924880510838743672432351;
1093
1094        rkT(1,2) = 1;
1095        rkT(1,2) = 1;
1096        rkT(1,3) = 0;
1097        rkT(2,1) = .4554100411010284672111720348287483;
1098        rkT(2,2) = -.6027050205505142336055860174143743;
1099        rkT(2,3) = -.4309321229203225731070721341350346;
1100        rkT(3,1) = 2.195823345445647152832799205549709;
1101        rkT(3,2) = -1.097911672722823576416399602774855;
1102        rkT(3,3) = .7850032632435902184104551358922130;
1103
1104        rkTinv(1,1) = .4205559181381766909344950150991349;
1105        rkTinv(1,2) = .3488903392193734304046467270632057;
1106        rkTinv(1,3) = .1915253879645878102698098373933487;
1107        rkTinv(2,1) = .5794440818618233090655049849008650;
1108        rkTinv(2,2) = -.3488903392193734304046467270632057;
1109        rkTinv(2,3) = -.1915253879645878102698098373933487;
1110        rkTinv(3,1) = -.3659705575742745254721332009249516;
1111        rkTinv(3,2) = -1.463882230297098101888532803699806;
1112        rkTinv(3,3) = .4702733607340189781407813565524989;
1113
1114        rkTinvAinv(1,1) = 1.104302803159744452668648155627548;
1115        rkTinvAinv(1,2) = .916122120694355522658740710823143;
1116        rkTinvAinv(1,3) = .5029105849749601702795812241441172;
1117        rkTinvAinv(2,1) = 1.895697196840255547331351844372453;
1118        rkTinvAinv(2,2) = 3.083877879305644477341259289176857;
1119        rkTinvAinv(2,3) = -1.502910584974960170279581224144117;
1120        rkTinvAinv(3,1) = .8362439183082935036129145574774502;
1121        rkTinvAinv(3,2) = -3.344975673233174014451658229909802;
1122        rkTinvAinv(3,3) = .312908409479233358005944466882642;
1123
1124        rkAinvT(1,1) = 2.625816818958466716011888933765282;
1125        rkAinvT(1,2) = 1.687091590520766641994055533117358;
1126        rkAinvT(1,3) = -2.508731754924880510838743672432351;
1127        rkAinvT(2,1) = 1.195823345445647152832799205549710;
1128        rkAinvT(2,2) = -2.097911672722823576416399602774855;
1129        rkAinvT(2,3) = .7850032632435902184104551358922130;
1130        rkAinvT(3,1) = 5.765829871932827589653709477334136;
1131        rkAinvT(3,2) = .1170850640335862051731452613329320;
1132        rkAinvT(3,3) = 4.078738281412060947659653944216779;
1133end;
1134return;
1135
1136        % /* Lobatto3C_Coefficients */
1137function Gauss_Coefficients()
1138       
1139        global rkMethod
1140        global rkT rkTinv rkTinvAinv rkAinvT rkA rkB rkC rkD rkE
1141        global rkBgam rkBhat rkTheta rkGamma rkAlpha rkBeta rkELO
1142        global GAU
1143
1144
1145   rkMethod = GAU;
1146   
1147   b0 = 0.1;
1148
1149   %/* The coefficients of the Gauss method */
1150   rkA(1,1) = .1388888888888888888888888888888889;
1151   rkA(1,2) = -.359766675249389034563954710966045e-01;
1152   rkA(1,3) = .97894440153083260495800422294756e-02;
1153   rkA(2,1) = .3002631949808645924380249472131556;
1154   rkA(2,2) = .2222222222222222222222222222222222;
1155   rkA(2,3) = -.224854172030868146602471694353778e-01;
1156   rkA(3,1) = .2679883337624694517281977355483022;
1157   rkA(3,2) = .4804211119693833479008399155410489;
1158   rkA(3,3) = .1388888888888888888888888888888889;
1159
1160   rkB(1) = .2777777777777777777777777777777778;
1161   rkB(2) = .4444444444444444444444444444444444;
1162   rkB(3) = .2777777777777777777777777777777778;
1163
1164   rkC(1) = .1127016653792583114820734600217600;
1165   rkC(2) = .5000000000000000000000000000000000;
1166   rkC(3) = .8872983346207416885179265399782400;
1167
1168%   /* Classical error estimator, embedded solution: */
1169   rkBhat(1) = b0;
1170   rkBhat(2) = -1.4788305577012361475298775666303999*b0...
1171                 +.27777777777777777777777777777777778;
1172   rkBhat(3) = .44444444444444444444444444444444444...
1173                 +.66666666666666666666666666666666667*b0;
1174   rkBhat(4) = -.18783610896543051913678910003626672*b0...
1175                 +.27777777777777777777777777777777778;
1176
1177 %  /* New solution: h Sum_j b_j f(Z_j) = sum d_j Z_j */
1178   rkD(1) = (.1666666666666666666666666666666667e+01);
1179   rkD(2) = (-.1333333333333333333333333333333333e+01);
1180   rkD(3) = (.1666666666666666666666666666666667e+01);
1181
1182   %/* Classical error estimator: */
1183   %/* H* Sum (B_j-Bhat_j)*f(Z_j) = H*E(0)*f(0) + Sum E_j*Z_j */
1184   rkE(1) = .2153144231161121782447335303806954*b0;
1185   rkE(2) = (-2.825278112319014084275808340593191)*b0;
1186   rkE(3) = .2870858974881495709929780405075939*b0;
1187   rkE(4) = (-.4558086256248162565397206448274867e-01)*b0;
1188
1189   %/* Sdirk error estimator */
1190   rkBgam(1) = 0;
1191   rkBgam(2) = .2373339543355109188382583162660537;
1192   rkBgam(3) = .5879873931885192299409334646982414;
1193   rkBgam(4) = -.4063577064014232702392531134499046e-01;
1194   rkBgam(5) = .2153144231161121782447335303806955;
1195
1196   %/* H* Sum Bgam_j*f(Z_j) = H*Bgam(0)*f(0) + Sum Theta_j*Z_j */
1197   rkTheta(1) = (-2.594040933093095272574031876464493);
1198   rkTheta(2) = 1.824611539036311947589425112250199;
1199   rkTheta(3) = .1856563166634371860478043996459493;
1200
1201   %/* ELO = local order of classical error estimator */
1202   rkELO = 4.0;
1203
1204   %/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1205   %~~~> Diagonalize the RK matrix:
1206   % rkTinv * inv(rkA) * rkT =
1207   %           |  rkGamma      0           0     |
1208   %           |      0      rkAlpha   -rkBeta   |
1209   %           |      0      rkBeta     rkAlpha  |
1210   %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1211
1212   rkGamma = 4.644370709252171185822941421408064;
1213   rkAlpha = 3.677814645373914407088529289295970;
1214   rkBeta  = 3.508761919567443321903661209182446;
1215
1216   rkT(1,1) = .7215185205520017032081769924397664e-01;
1217   rkT(1,2) = -.8224123057363067064866206597516454e-01;
1218   rkT(1,3) = -.6012073861930850173085948921439054e-01;
1219   rkT(2,1) = .1188325787412778070708888193730294;
1220   rkT(2,2) = .5306509074206139504614411373957448e-01;
1221   rkT(2,3) = .3162050511322915732224862926182701;
1222   rkT(3,1) = 1;
1223   rkT(3,2) = 1;
1224   rkT(3,3) = 0;
1225
1226   rkTinv(1,1) = 5.991698084937800775649580743981285;
1227   rkTinv(1,2) = 1.139214295155735444567002236934009;
1228   rkTinv(1,3) = .4323121137838583855696375901180497;
1229   rkTinv(2,1) = -5.991698084937800775649580743981285;
1230   rkTinv(2,2) = -1.139214295155735444567002236934009;
1231   rkTinv(2,3) = .5676878862161416144303624098819503;
1232   rkTinv(3,1) = -1.246213273586231410815571640493082;
1233   rkTinv(3,2) = 2.925559646192313662599230367054972;
1234   rkTinv(3,3) = -.2577352012734324923468722836888244;
1235
1236   rkTinvAinv(1,1) = 27.82766708436744962047620566703329;
1237   rkTinvAinv(1,2) = 5.290933503982655311815946575100597;
1238   rkTinvAinv(1,3) = 2.007817718512643701322151051660114;
1239   rkTinvAinv(2,1) = (-17.66368928942422710690385180065675);
1240   rkTinvAinv(2,2) = (-14.45491129892587782538830044147713);
1241   rkTinvAinv(2,3) = 2.992182281487356298677848948339886;
1242   rkTinvAinv(3,1) = (-25.60678350282974256072419392007303);
1243   rkTinvAinv(3,2) = 6.762434375611708328910623303779923;
1244   rkTinvAinv(3,3) = 1.043979339483109825041215970036771;
1245
1246   rkAinvT(1,1) = .3350999483034677402618981153470483;
1247   rkAinvT(1,2) = -.5134173605009692329246186488441294;
1248   rkAinvT(1,3) = .6745196507033116204327635673208923e-01;
1249   rkAinvT(2,1) = .5519025480108928886873752035738885;
1250   rkAinvT(2,2) = 1.304651810077110066076640761092008;
1251   rkAinvT(2,3) = .9767507983414134987545585703726984;
1252   rkAinvT(3,1) = 4.644370709252171185822941421408064;
1253   rkAinvT(3,2) = 3.677814645373914407088529289295970;
1254   rkAinvT(3,3) = -3.508761919567443321903661209182446;
1255return;
1256 %/* Gauss_Coefficients */
1257
1258%/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1259%       The coefficients of the 3-stage Gauss method
1260%       (given to ~30 accurate digits)
1261%       The parameter b3 can be chosen by the user
1262%       to tune the error estimator
1263%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1264function Radau1A_Coefficients()
1265        global rkMethod
1266        global rkT rkTinv rkTinvAinv rkAinvT rkA rkB rkC rkD rkE
1267        global rkBgam rkBhat rkTheta rkGamma rkAlpha rkBeta rkELO
1268        global  R1A
1269
1270        global ZERO
1271
1272
1273%   /* The coefficients of the Radau1A method */
1274   b0 = 0.1;
1275   rkMethod = R1A;
1276   rkA(1,1) =  .1111111111111111111111111111111111;
1277   rkA(1,2) = (-.1916383190435098943442935597058829);
1278   rkA(1,3) = (.8052720793239878323318244859477174e-01);
1279   rkA(2,1) = .1111111111111111111111111111111111;
1280   rkA(2,2) = .2920734116652284630205027458970589;
1281   rkA(2,3) = (-.481334970546573839513422644787591e-01);
1282   rkA(3,1) = .1111111111111111111111111111111111;
1283   rkA(3,2) = .5370223859435462728402311533676479;
1284   rkA(3,3) = .1968154772236604258683861429918299;
1285
1286   rkB(1) = .1111111111111111111111111111111111;
1287   rkB(2) = .5124858261884216138388134465196080;
1288   rkB(3) = .3764030627004672750500754423692808;
1289
1290   rkC(1) = 0;
1291   rkC(2) = .3550510257216821901802715925294109;
1292   rkC(3) = .8449489742783178098197284074705891;
1293
1294  % /* Classical error estimator, embedded solution: */
1295   rkBhat(1) = b0;
1296   rkBhat(2) = .11111111111111111111111111111111111-b0;
1297   rkBhat(3) = .51248582618842161383881344651960810;
1298   rkBhat(4) = .37640306270046727505007544236928079;
1299
1300   %/* New solution: H* Sum B_j*f(Z_j) = Sum D_j*Z_j */
1301   rkD(1) = .3333333333333333333333333333333333;
1302   rkD(2) = -.8914115380582557157653087040196127;
1303   rkD(3) = 1.558078204724922382431975370686279;
1304
1305   %/* Classical error estimator: */
1306   %/* H* Sum (b_j-bhat_j) f(Z_j) = H*E(0)*F(0) + Sum E_j Z_j */
1307   rkE(1) = .2748888295956773677478286035994148*b0;
1308   rkE(2) = -1.374444147978386838739143017997074*b0;
1309   rkE(3) = -1.335337922441686804550326197041126*b0;
1310   rkE(4) = .235782604058977333559011782643466*b0;
1311
1312   %/* Sdirk error estimator */
1313   rkBgam(1) = ZERO;
1314   rkBgam(2) = (.1948150124588532186183490991130616e-01);
1315   rkBgam(3) = .7575249005733381398986810981093584;
1316   rkBgam(4) = (-.518952314149008295083446116200793e-01);
1317   rkBgam(5) = .2748888295956773677478286035994148;
1318
1319   %/* H* Sum Bgam_j*f(Z_j) = H*Bgam(0)*f(0) + Sum Theta_j*Z_j */
1320   rkTheta(1) = (-1.224370034375505083904362087063351);
1321   rkTheta(2) = .9340045331532641409047527962010133;
1322   rkTheta(3) = .4656990124352088397561234800640929;
1323
1324   %/* ELO = local order of classical error estimator */
1325   rkELO = 4.0;
1326
1327   %/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1328   %~~~> Diagonalize the RK matrix:
1329   % rkTinv * inv(rkA) * rkT =
1330   %           |  rkGamma      0           0     |
1331   %           |      0      rkAlpha   -rkBeta   |
1332   %           |      0      rkBeta     rkAlpha  |
1333   %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1334
1335   rkGamma = 3.637834252744495732208418513577775;
1336   rkAlpha = 2.681082873627752133895790743211112;
1337   rkBeta  = 3.050430199247410569426377624787569;
1338
1339   rkT(1,1) = .424293819848497965354371036408369;
1340   rkT(1,2) = -.3235571519651980681202894497035503;
1341   rkT(1,3) = -.522137786846287839586599927945048;
1342   rkT(2,1) = .57594609499806128896291585429339e-01;
1343   rkT(2,2) = .3148663231849760131614374283783e-02;
1344   rkT(2,3) = .452429247674359778577728510381731;
1345   rkT(3,1) = 1;
1346   rkT(3,2) = 1;
1347   rkT(3,3) = 0;
1348
1349   rkTinv(1,1) = 1.233523612685027760114769983066164;
1350   rkTinv(1,2) = 1.423580134265707095505388133369554;
1351   rkTinv(1,3) = .3946330125758354736049045150429624;
1352   rkTinv(2,1) = -1.233523612685027760114769983066164;
1353   rkTinv(2,2) = -1.423580134265707095505388133369554;
1354   rkTinv(2,3) = .6053669874241645263950954849570376;
1355   rkTinv(3,1) = -.1484438963257383124456490049673414;
1356   rkTinv(3,2) = 2.038974794939896109682070471785315;
1357   rkTinv(3,3) = -.544501292892686735299355831692542e-01;
1358
1359   rkTinvAinv(1,1) = 4.487354449794728738538663081025420;
1360   rkTinvAinv(1,2) = 5.178748573958397475446442544234494;
1361   rkTinvAinv(1,3) = 1.435609490412123627047824222335563;
1362   rkTinvAinv(2,1) = -2.854361287939276673073807031221493;
1363   rkTinvAinv(2,2) = -1.003648660720543859000994063139137e+01;
1364   rkTinvAinv(2,3) = 1.789135380979465422050817815017383;
1365   rkTinvAinv(3,1) = -4.160768067752685525282947313530352;
1366   rkTinvAinv(3,2) = 1.124128569859216916690209918405860;
1367   rkTinvAinv(3,3) = 1.700644430961823796581896350418417;
1368
1369   rkAinvT(1,1) = 1.543510591072668287198054583233180;
1370   rkAinvT(1,2) = -2.460228411937788329157493833295004;
1371   rkAinvT(1,3) = -.412906170450356277003910443520499;
1372   rkAinvT(2,1) = .209519643211838264029272585946993;
1373   rkAinvT(2,2) = 1.388545667194387164417459732995766;
1374   rkAinvT(2,3) = 1.20339553005832004974976023130002;
1375   rkAinvT(3,1) = 3.637834252744495732208418513577775;
1376   rkAinvT(3,2) = 2.681082873627752133895790743211112;
1377   rkAinvT(3,3) = -3.050430199247410569426377624787569;
1378return;
1379% /* Radau1A_Coefficients */
1380
1381%/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1382%       The coefficients of the 4-stage Lobatto-3A method
1383%       (given to ~30 accurate digits)
1384%  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1385function Lobatto3A_Coefficients()
1386       
1387        global rkMethod
1388        global rkT rkTinv rkTinvAinv rkAinvT rkA rkB rkC rkD rkE
1389        global rkBgam rkBhat rkTheta rkF rkGamma rkAlpha rkBeta rkELO
1390        global L3A
1391
1392
1393   %/* The coefficients of the Lobatto-3A method */
1394   rkMethod = L3A;
1395
1396   rkA(1,1) = 0;
1397   rkA(1,2) = 0;
1398   rkA(1,3) = 0;
1399   rkA(1,4) = 0;
1400   rkA(2,1) = .11030056647916491413674311390609397;
1401   rkA(2,2) = .1896994335208350858632568860939060;
1402   rkA(2,3) = -.339073642291438837776604807792215e-01;
1403   rkA(2,4) = .1030056647916491413674311390609397e-01;
1404   rkA(3,1) = .73032766854168419196590219427239365e-01;
1405   rkA(3,2) = .4505740308958105504443271474458881;
1406   rkA(3,3) = .2269672331458315808034097805727606;
1407   rkA(3,4) = -.2696723314583158080340978057276063e-01;
1408   rkA(4,1) = .83333333333333333333333333333333333e-01;
1409   rkA(4,2) = .4166666666666666666666666666666667;
1410   rkA(4,3) = .4166666666666666666666666666666667;
1411   rkA(4,4) = .8333333333333333333333333333333333e-01;
1412
1413   rkB(1) = .83333333333333333333333333333333333e-01;
1414   rkB(2) = .4166666666666666666666666666666667;
1415   rkB(3) = .4166666666666666666666666666666667;
1416   rkB(4) = .8333333333333333333333333333333333e-01;
1417
1418   rkC(1) = 0;
1419   rkC(2) = .2763932022500210303590826331268724;
1420   rkC(3) = .7236067977499789696409173668731276;
1421   rkC(4) = 1;
1422
1423   %/* New solution: H*Sum B_j*f(Z_j) = Sum D_j*Z_j */
1424   rkD(1) = 0;
1425   rkD(2) = 0;
1426   rkD(3) = 0;
1427   rkD(4) = 1;
1428
1429   %/* Classical error estimator, embedded solution: */
1430   rkBhat(1) = .90909090909090909090909090909090909e-01;
1431   rkBhat(2) = .39972675774621371442114262372173276;
1432   rkBhat(3) = .43360657558711961891219070961160058;
1433   rkBhat(4) = .15151515151515151515151515151515152e-01;
1434
1435   %/* Classical error estimator: */
1436   %/* H* Sum (B_j-Bhat_j)*f(Z_j) = H*E(0)*f(0) + Sum E_j*Z_j */
1437
1438   rkE(1) = .1957403846510110711315759367097231e-01;
1439   rkE(2) = -.1986820345632580910316020806676438;
1440   rkE(3) = .1660586371214229125096727578826900;
1441   rkE(4) = -.9787019232550553556578796835486154e-01;
1442
1443   %/* Sdirk error estimator: */
1444   rkF(1) = 0;
1445   rkF(2) = -.66535815876916686607437314126436349;
1446   rkF(3) = 1.7419302743497277572980407931678409;
1447   rkF(4) = -1.2918865386966730694684011822841728;
1448
1449   %/* ELO = local order of classical error estimator */
1450   rkELO = 4.0;
1451
1452   %/* Sdirk error estimator: */
1453   rkBgam(1) = .2950472755430528877214995073815946e-01;
1454   rkBgam(2) = .5370310883226113978352873633882769;
1455   rkBgam(3) = .2963022450107219354980459699450564;
1456   rkBgam(4) = -.7815248400375080035021681445218837e-01;
1457   rkBgam(5) = .2153144231161121782447335303806956;
1458
1459   %/* H* Sum Bgam_j*f(Z_j) = H*Bgam(0)*f(0) + Sum Theta_j*Z_j */
1460   rkTheta(1) = 0.0;
1461   rkTheta(2) = -.6653581587691668660743731412643631;
1462   rkTheta(3) = 1.741930274349727757298040793167842;
1463   rkTheta(4) = -.291886538696673069468401182284174;
1464
1465
1466   %/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1467   %~~~> Diagonalize the RK matrix:
1468   % rkTinv * inv(rkA) * rkT =
1469   %          |  rkGamma      0           0     |
1470   %           |      0      rkAlpha   -rkBeta   |
1471   %           |      0      rkBeta     rkAlpha  |
1472   %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1473
1474   rkGamma = 4.644370709252171185822941421408063;
1475   rkAlpha = 3.677814645373914407088529289295968;
1476   rkBeta  = 3.508761919567443321903661209182446;
1477
1478   rkT(1,1) = .5303036326129938105898786144870856e-01;
1479   rkT(1,2) = -.7776129960563076320631956091016914e-01;
1480   rkT(1,3) = .6043307469475508514468017399717112e-02;
1481   rkT(2,1) = .2637242522173698467283726114649606;
1482   rkT(2,2) = .2193839918662961493126393244533346;
1483   rkT(2,3) = .3198765142300936188514264752235344;
1484   rkT(3,1) = 1;
1485   rkT(3,2) = 0;
1486   rkT(3,3) = 0;
1487
1488   rkTinv(1,1) = 7.695032983257654470769069079238553;
1489   rkTinv(1,2) = -.1453793830957233720334601186354032;
1490   rkTinv(1,3) = .6302696746849084900422461036874826;
1491   rkTinv(2,1) = -7.695032983257654470769069079238553;
1492   rkTinv(2,2) = .1453793830957233720334601186354032;
1493   rkTinv(2,3) = .3697303253150915099577538963125174;
1494   rkTinv(3,1) = -1.066660885401270392058552736086173;
1495   rkTinv(3,2) = 3.146358406832537460764521760668932;
1496   rkTinv(3,3) = -.7732056038202974770406168510664738;
1497
1498   rkTinvAinv(1,1) = 35.73858579417120341641749040405149;
1499   rkTinvAinv(1,2) = -.675195748578927863668368190236025;
1500   rkTinvAinv(1,3) = 2.927206016036483646751158874041632;
1501   rkTinvAinv(2,1) = -24.55824590667225493437162206039511;
1502   rkTinvAinv(2,2) = -10.50514413892002061837750015342036;
1503   rkTinvAinv(2,3) = 4.072793983963516353248841125958369;
1504   rkTinvAinv(3,1) = -30.92301972744621647251975054630589;
1505   rkTinvAinv(3,2) = 12.08182467154052413351908559269928;
1506   rkTinvAinv(3,3) = -1.546411207640594954081233702132946;
1507
1508   rkAinvT(1,1) = .2462926658317812882584158369803835;
1509   rkAinvT(1,2) = -.2647871194157644619747121197289574;
1510   rkAinvT(1,3) = .2950720515900466654896406799284586;
1511   rkAinvT(2,1) = 1.224833192317784474576995878738004;
1512   rkAinvT(2,2) = 1.929224190340981580557006261869763;
1513   rkAinvT(2,3) = .4066803323234419988910915619080306;
1514   rkAinvT(3,1) = 4.644370709252171185822941421408064;
1515   rkAinvT(3,2) = 3.677814645373914407088529289295968;
1516   rkAinvT(3,3) = -3.508761919567443321903661209182446;
1517return;
1518                %/* Lobatto3A_Coefficients */
1519
Note: See TracBrowser for help on using the repository browser.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |