source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int/sdirk.m @ 3961

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

Merge of branch palm4u into trunk

File size: 28.0 KB
Line 
1function [T,Y,RCNTRL,ICNTRL,RSTATUS,ISTATUS]=...
2    SdirkInt(Function,Tspan,Y0,Options,RCNTRL,ICNTRL)
3%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
4%  SDIRK - Implementation of several SDIRK methods:
5%               * SDIRK4a                                               
6%               * SDIRK4b                                               
7%               * SDIRK3a                                               
8%               * SDIRK2b                                             
9%               * SDIRK2a
10%                                                                     
11%    Solves the system y'=F(t,y) using a SDIRK Method
12%                                                                     
13%
14%    For details on SDIRK 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%      RSTATUS     - real output parameters (explained below)
52%      ISTATUS     - 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%
62%    ICNTRL(1) = 0: AbsTol, RelTol are N-dimensional vectors
63%              = 1: AbsTol, RelTol are scalars
64%
65%    ICNTRL(2)  -> selection of a particular SDIRK method
66%        = 0 :    Sdirk2a(Default)
67%        = 1 :    Sdirk2a
68%        = 2 :    Sdirk2b
69%        = 3 :    Sdirk3a
70%        = 4 :    Sdirk4a
71%        = 5 :    Sdirk4b
72%       
73%       
74%
75%    ICNTRL(3)  -> maximum number of integration steps
76%        For ICNTRL(3)=0) the default value of 100000 is used
77%
78%    RCNTRL(1)  -> Hmin, lower bound for the integration step size
79%          It is strongly recommended to keep Hmin = ZERO
80%    RCNTRL(2)  -> Hmax, upper bound for the integration step size
81%    RCNTRL(3)  -> Hstart, starting value for the integration step size
82%
83%    RCNTRL(4)  -> FacMin, lower bound on step decrease factor (default=0.2)
84%    RCNTRL(5)  -> FacMax, upper bound on step increase factor (default=6)
85%    RCNTRL(6)  -> FacRej, step decrease factor after multiple rejections
86%                          (default=0.1)
87%    RCNTRL(7)  -> FacSafe, by which the new step is slightly smaller
88%         than the predicted value  (default=0.9)
89%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
90%
91%    RSTAT and ISTAT on output:
92%
93%    ISTATUS(1)  -> No. of function calls
94%    ISTATUS(2)  -> No. of jacobian calls
95%    ISTATUS(3)  -> No. of steps
96%    ISTATUS(4)  -> No. of accepted steps
97%    ISTATUS(5)  -> No. of rejected steps (except at very beginning)
98%    ISTATUS(6)  -> No. of LU decompositions
99%    ISTATUS(7)  -> No. of forward/backward substitutions
100%    ISTATUS(8)  -> No. of singular matrix decompositions
101%
102%    RSTATUS(1)  -> Texit, the time corresponding to the
103%                     computed Y upon return
104%    RSTATUS(2)  -> Hexit, last accepted step before exit
105%    RSTATUS(3)  -> Hnew, last predicted step (not yet taken)
106%                   For multiple restarts, use Hnew as Hstart
107%                     in the subsequent run
108%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
109global Nfun Njac Nstp Nacc Nrej Ndec Nsol Nsng Ntexit Nhexit Nhnew
110% Parse ODE options
111AbsTol = odeget(Options, 'AbsTol');
112if isempty(AbsTol)
113    AbsTol = 1.0e-3;
114end
115Hstart = odeget(Options, 'InitialStep');
116if isempty(Hstart)
117    Hstart = 0;
118end
119Jacobian = odeget(Options, 'Jacobian');
120if isempty(Jacobian)
121    error('A Jacobian function is required.');
122end
123Hmax = odeget(Options, 'MaxStep');
124if isempty(Hmax)
125    Hmax = 0;
126end
127RelTol = odeget(Options, 'RelTol');
128if isempty(RelTol)
129    RelTol = 1.0e-4;
130end
131% Initialize statistics
132Nfun=1; Njac=2; Nstp=3; Nacc=4; Nrej=5; Ndec=6; Nsol=7; Nsng=8;
133Ntexit=1; Nhexit=2; Nhnew=3;
134RSTATUS = zeros(20,1);
135ISTATUS = zeros(20,1);
136
137% Get problem size
138steps = length(Tspan);
139N = max(size(Y0));
140
141% Initialize tolerances
142ATOL = ones(N,1)*AbsTol;
143RTOL = ones(N,1)*RelTol;
144
145% Initialize step
146RCNTRL(2) = max(0, Hmax);
147RCNTRL(3) = max(0, Hstart);
148
149% Integrate
150Y = zeros(N,steps);
151T = zeros(steps,1);
152Y(:,1) = Y0;
153T(1) = Tspan(1);
154t=cputime;
155for i=2:steps
156    [T(i), Y(:,i), IERR,ISTATUS,RSTATUS] = ...
157        Sdirk_Integrate(N, Y(:,i-1), T(i-1), Tspan(i), ...
158            Function, Jacobian, ATOL, RTOL, RCNTRL, ICNTRL,ISTATUS,RSTATUS);
159
160    if IERR < 0
161        error(['SDIRK exited with IERR=',num2str(IERR)]);
162    end
163
164end
165cputime-t
166
167
168
169function [T,Y,Ierr,ISTATUS,RSTATUS] = Sdirk_Integrate(N,Y0,TIN, TOUT,OdeFunction1,...
170    OdeJacobian1,ATOL,RTOL,RCNTRL,ICNTRL,ISTATUS,RSTATUS)
171global ZERO ONE Nfun Njac Nstp Nacc Nrej Ndec Nsol Nsng Ntexit Nhexit Nhnew
172global Smax S2A S2B S3A S4A S4B
173global rkTheta rkAlpha
174ZERO = 0; ONE=1.0; Nfun = 1; Njac = 2; Nstp = 3; Nacc = 4; Nrej = 5;Ndec=6;
175Nsol = 7; Nsng = 8; Ntexit = 1; Nhexit = 2; Nhnew = 3; Smax = 5; S2A = 1;
176S2B = 2; S3A = 3; S4A = 4; S4B = 5;
177rkTheta = zeros(Smax);
178rkAlpha = zeros(Smax);
179
180%Local Variables
181ISTATUS = zeros(20,1);
182RSTATUS = zeros(20,1);
183
184% Fine-Tune the integrator
185T1 = TIN;
186T2 = TOUT;
187RTOL=ones(size(Y0));
188RTOL=1e-3*RTOL;
189ATOL=RTOL;
190
191    [T,Y,Ierr,RSTATUS,ISTATUS] = SDIRK(N, T1,T2,Y0,RTOL,ATOL,RCNTRL,ICNTRL,RSTATUS,ISTATUS,OdeFunction1,OdeJacobian1);
192
193if Ierr < 0
194    disp(['SDIRK:Unsuccessful exit at ',num2str(TIN), ' ',num2str(Ierr)]);
195end
196return;
197
198
199
200
201function [T, Y, Ierr,RSTATUS,ISTATUS] = SDIRK(N,T1,T2,Y,AbsTol,RelTol,...
202                        RCNTRL,ICNTRL,RSTATUS,ISTATUS,OdeFunction1,OdeJacobian1)
203
204global ZERO ONE
205
206
207Max_no_steps = 0;Hmin = 0;Hmax = 0; Hstart = 0; Roundoff = eps; FacMin = 0;
208FacMax = 0; FacSafe = 0; FacRej = 0;
209if(ICNTRL(1)==0)
210    ITOL = 1;
211else
212    ITOL = 0;
213end
214Ierr = 0;
215%~~~> Method Selection
216switch (ICNTRL(2))
217   
218    case (0)
219        Sdirk2a();
220    case (1)
221        Sdirk2a();
222    case (2)
223        Sdirk2b();
224    case(3)
225        Sdirk3a();
226    case(4)
227        Sdirk4a();
228    case(5)
229        Sdirk4b();
230    otherwise
231        Sdirk2a();
232end
233
234if (ICNTRL(3)==ZERO)
235    Max_no_steps = 200000;
236end
237if(ICNTRL(3) > ZERO)
238        Max_no_steps = ICNTRL(3);
239end
240
241if(ICNTRL(3) < ZERO)
242    disp(['User-Selected ICNTRL(3)=',num2str(ICNTRL(3))]);
243    SDIRK_ErrorMsg(-1,T1,ZERO,Ierr);
244end;
245%~~>StartNewton Extrapolate for starting values of Newton iterations
246if(ICNTRL(5)==ZERO)
247    StartNewton = 1;
248else
249    StartNewton = 0;
250end;
251
252if(ICNTRL(4)>=ZERO)
253    if(ICNTRL(4)==ZERO)
254        NewtonMaxit = 8;
255    else
256        NewtonMaxit = ICNTRL(4);
257    end
258else
259    disp(['User-Selected NewtonMaxit ICNTRL(4)=',num2str(ICNTRL(4))]);
260    SDIRK_ErrorMsg(-2,T1,ZERO,Ierr);
261end;
262
263%~~> Lower bound on the step size (Positive value)
264
265if(RCNTRL(1)==ZERO)
266    Hmin = ZERO;
267end
268if (RCNTRL(1)>ZERO)
269        Hmin = RCNTRL(1);
270end   
271if(RCNTRL(1)<ZERO)
272    disp(['User-Selected Hmin, RCNTRL(1)=',num2str(RCNTRL(1))]);
273    SDIRK_ErrorMsg(-3,T1,ZERO,Ierr);
274end;
275% Upper bound on the step size : (Positive value)
276
277if(RCNTRL(2)==ZERO)
278    Hmax = abs(T2-T1);
279end
280if(RCNTRL(2)>ZERO)
281    Hmax = min(abs(T2-T1),RCNTRL(2));
282end;
283if(RCNTRL(2)<ZERO)
284    disp(['User-Selected Hmax, RCNTRL(2)=',num2str(RCNTRL(2))]);
285    SDIRK_ErrorMsg(-3,T1,ZERO,Ierr);
286end;
287%~~> Starting step size
288
289if(RCNTRL(3)>=ZERO)
290    if(RCNTRL(3)==ZERO)
291        Hstart = max(Hmin,Roundoff);
292    else
293        Hstart = min(abs(RCNTRL(3)),abs(T2-T1));
294    end
295else
296    disp(['User-Selected Hstart, RCNTRL(3)=',num2str(RCNTRL(3))]);
297    SDIRK_ErrorMsg(-3,T1,ZERO,Ierr);
298end;
299%~~> Stepsize can be changed such that FacMin < Hnew/Hexit < FacMax
300if(RCNTRL(4) >= ZERO)
301    if(RCNTRL(4) == ZERO)
302        FacMin = 0.2;
303    else
304        FacMin = RCNTRL(4);
305    end
306else
307    disp(['User-Selected FacMin, RCNTRL(4)=',num2str(RCNTRL(4))]);
308    SDIRK_ErrorMsg(-4,T1,ZERO,Ierr);
309end
310
311
312if (RCNTRL(5)>=ZERO)
313    if(RCNTRL(5)==ZERO)
314        FacMax = 10.0;
315    else
316        FacMax = RCNTRL(5);
317    end
318else
319    disp(['User-Selected FacMax, RCNTRL(5)=',num2str(RCNTRL(5))]);
320    SDIRK_ErrorMsg(-4,T1,ZERO,Ierr);
321end
322
323%~~> FacRej Factor to decrease step after 2 successive rejections
324
325if(RCNTRL(6)>=ZERO)
326    if(RCNTRL(6) == ZERO)
327        FacRej = 0.1;
328    else
329        FacRej = RCNTRL(6);
330    end
331else
332    disp(['User-Selected FacRej, RCNTRL(6)=',num2str(RCNTRL(6))]);
333    SDIRK_ErrorMsg(-4,T1,ZERO,Ierr);
334end;
335
336%~~> FacSafe Safety Factor in the computation of new step size
337
338if(RCNTRL(7)>=ZERO)
339    if(RCNTRL(7)==ZERO)
340        FacSafe = 0.9;
341    else
342        FacSafe = RCNTRL(7);
343    end
344else
345    disp(['User-Selected FacSafe, RCNTRL(7)=',num2str(RCNTRL(7))]);
346    SDIRK_ErrorMsg(-4,T1,ZERO,Ierr);
347end
348
349%~~> ThetaMin decides whether the Jacobian should be recomputed
350
351if(RCNTRL(8) == ZERO)
352    ThetaMin = 1.0E-3;
353else
354    ThetaMin = RCNTRL(8);
355end;
356
357%~~> Stopping Criteria for Newtons method
358
359if(RCNTRL(9) == ZERO)
360    NewtonTol = 3.0e-4;
361else
362    NewtonTol = RCNTRL(9);
363end;
364
365%~~> Qmin, Qmax : IF Qmin <Hnew/Hold<Qmax STEP SIZE = constant
366
367if(RCNTRL(10) == ZERO)
368    Qmin = ONE;
369else
370    Qmin = RCNTRL(10);
371end
372
373if(RCNTRL(11) == ZERO)
374    Qmax = 1.2;
375else
376    Qmax = RCNTRL(11);
377end
378
379%~~> Check if tolerances are reasonable
380if(ITOL == ZERO)
381    if((AbsTol(1)<=ZERO) || RelTol(1) <=10*Roundoff)
382        SDIRK_ErrorMsg(-5,T1,ZERO,Ierr);
383    end;
384else
385    for i=1:N
386        if((AbsTol(i)<=ZERO)||RelTol(i)<=10*Roundoff)
387            SDIRK_ErrorMsg(-5,T1,ZERO,Ierr);
388        end
389    end
390end
391if(Ierr < 0)
392    return;
393end
394
395
396
397[T,Y,Ierr,ISTATUS,RSTATUS] = SDIRK_Inetgrator(N,T1,T2,Y,Hstart,Hmin,...
398                            Hmax,Roundoff,AbsTol,RelTol,...
399                            ITOL,Max_no_steps,StartNewton,NewtonTol,ThetaMin,...
400                            FacSafe,FacMin,FacMax,FacRej,Qmin,Qmax,NewtonMaxit,ISTATUS,RSTATUS,Ierr,OdeFunction1,OdeJacobian1);
401 
402
403
404return;
405
406function[T,Y,Ierr,ISTATUS,RSTATUS] = SDIRK_Inetgrator(N,T1,T2,Y,Hstart,Hmin,...
407                            Hmax,Roundoff,AbsTol,RelTol,...
408                            ITOL,Max_no_steps,StartNewton,NewtonTol,ThetaMin,...
409                            FacSafe,FacMin,FacMax,FacRej,Qmin,Qmax,NewtonMaxit,ISTATUS,RSTATUS,Ierr,OdeFunction1,OdeJacobian1)
410global ZERO ONE Nfun Nstp Nacc Nrej Ntexit Nhexit Nhnew
411global rkS rkGamma rkELO
412global rkC rkD rkE rkTheta rkAlpha
413T =T1;
414Tdirection = sign(T2-T1);
415H = max(abs(Hmin),abs(Hstart));
416if(H<=10.0*Roundoff)
417    H = 1.0e-6;
418end;
419H = min(H,Hmax);
420H=sign(Tdirection)*H;
421SkipLU = 0;
422SkipJac = 0;
423Reject = 0;
424FirstStep = 1;
425CycleTloop = 0;
426SCAL = SDIRK_ErrorScale(N,ITOL,AbsTol,RelTol,Y);
427FJAC = zeros(N);
428while((T2-T)*Tdirection - Roundoff > ZERO)
429    if(SkipLU == 0)
430        [E1L,E1U,ISTATUS,ISING]=PrepareMatrix(N,H,T,Y,FJAC,SkipJac, SkipLU,Reject,ISTATUS,OdeJacobian1);
431        if(ISING ~= 0)
432            SDIRK_ErrorMsg(-8,T,H,Ierr);
433            return;
434        end;
435    end;
436    if(ISTATUS(Nstp)>Max_no_steps)
437        SDIRK_ErrorMsg(-6, T, H, Ierr);
438        return;
439    end;
440   
441    if(T+0.1*H==T || abs(H)<=Roundoff)
442        SDIRK_ErrorMsg(-7,T,H,Ierr);
443        return;
444    end
445    %~~> Simplified Newton Iterations
446    Z=zeros(rkS,N);
447    for istage=1:rkS
448        G=zeros(N,1);
449        if(istage>1)
450            aa=rkTheta(1:istage,istage);
451            G=Z(1:istage,:)'*aa;
452            if(StartNewton==1)
453                for j=1:istage
454                    Z(istage,:)=Z(istage,:)+rkAlpha(j,istage)*Z(j,:);
455                end
456               
457            end
458               
459                end
460                NewtonDone = 0;
461                Fac = 0.5;
462                for     NewtonIter=0:NewtonMaxit-1
463                        TMP=Z(istage,:)'+Y;
464                        RHS=OdeFunction1(T+rkC(istage)*H,TMP);
465                        ISTATUS(Nfun)=ISTATUS(Nfun)+1;
466                        RHS=H*rkGamma*RHS;
467                        RHS=RHS-Z(istage,:)';
468            RHS=RHS+G;
469                        [RHS,ISTATUS]=SDIRK_Solve(H,N,E1L,E1U,RHS,ISTATUS);
470                        NewtonIncrement = SDIRK_ErrorNorm(N,RHS,SCAL);
471                        if(NewtonIter==0)
472                                Theta = abs(ThetaMin);
473                                NewtonRate = 2.0;
474                        else
475                                Theta = NewtonIncrement/NewtonIncrementOld;
476                                if (Theta < 0.99)
477                                        NewtonRate = Theta/(1-Theta);
478                                        NewtonPredictedErr = NewtonIncrement*Theta^((NewtonMaxit - (NewtonIter+1))/(ONE-Theta));
479                                        if(NewtonPredictedErr>=NewtonTol)
480                                                Qnewton = min(10,NewtonPredictedErr/NewtonTol);
481                                                Fac = 0.8*Qnewton ^(-ONE/(NewtonMaxit-NewtonIter));
482                                                break;
483                                        end
484                                else
485                                        break;
486                                end
487                        end;
488                        NewtonIncrementOld = NewtonIncrement;
489                        Z(istage,:)=Z(istage,:)+RHS';
490                        NewtonDone=(NewtonRate*NewtonIncrement<=NewtonTol);
491                        if(NewtonDone == 1)
492                                break;
493                        end;
494                end;
495
496                if(NewtonDone == 0)
497                        H=Fac*H;
498                        Reject =1;
499                        SkipJac = 1;
500                        SkipLU = 0;
501                        CycleTloop = 1;
502                end;
503                if(CycleTloop == 1)
504                        CycleTloop = 0;
505                        break;
506                end
507        end;
508        if(CycleTloop==0)
509        %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
510        %~~~~ Error Estimation ~~~~~~~~~~~~~~%
511        %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
512        ISTATUS(Nstp) = ISTATUS(Nstp)+1;
513        TMP = zeros(N,1);
514    rr=rkE(1:rkS);
515    TMP = (rr*Z)';
516       
517        [TMP,ISTATUS]=  SDIRK_Solve(H,N,E1L,E1U,TMP,ISTATUS);
518        Err = SDIRK_ErrorNorm (N,TMP,SCAL);
519        %/*~~~~> Computation of new step size Hnew */
520
521        Fac = FacSafe *Err^(-ONE/rkELO);
522        Fac = max (FacMin, min(FacMax,Fac));
523    Hnew = H*Fac;
524        %/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
525        %/*~~~>  Accept/Reject step    */
526        %/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
527
528        if(Err < ONE)
529                FirstStep = 0;
530                ISTATUS(Nacc)=ISTATUS(Nacc)+1;
531                T=T+H;
532        Y=Y+(rkD(1:rkS)*Z)';
533                SCAL = SDIRK_ErrorScale(N,ITOL,AbsTol,RelTol,Y);
534                Hnew = Tdirection*min(abs(Hnew),Hmax);
535                RSTATUS(Ntexit)=T;
536                RSTATUS(Nhexit)=H;
537                RSTATUS(Nhnew)=Hnew;
538                if(Reject == 1)
539                        Hnew = Tdirection*min(abs(Hnew),abs(H));
540                end;
541                Reject = 0;
542                if((T+Hnew/Qmin - T2)*Tdirection > ZERO)
543                        H=T2-T;
544                else
545                        Hratio = Hnew/H;
546                        SkipLU = ((Theta<=ThetaMin) && (Hratio >=Qmin) && (Hratio<=Qmax));
547                        if(SkipLU==0)
548                                H=Hnew;
549                        end
550                end
551                SkipJac = (Theta<=ThetaMin);
552                SkipJac = 0;
553        else
554                if(FirstStep==1 || Reject==1)
555                        H=FacRej*H;
556                else
557                        H=Hnew;
558                end;
559                Reject=1;
560                SkipJac =1;
561                SkipLU=0;
562                if(ISTATUS(Nacc)>=1)
563                        ISTATUS(Nrej)=ISTATUS(Nrej)+1;
564                end
565        end
566        end
567end
568
569Ierr =1;
570return;
571
572
573
574
575function [SCAL] = SDIRK_ErrorScale(N,ITOL,AbsTol,RelTol,Y)
576ZERO =0;
577if(ITOL==ZERO)
578    temp = RelTol(1)*Y;
579    temp =(AbsTol(1)+abs(temp));
580    SCAL=1./temp;
581else
582    SCAL=1./(AbsTol+RelTol.*abs(Y));
583end
584return;
585
586
587function [Err] = SDIRK_ErrorNorm(N,Y,SCAL)
588ZERO=0;
589temp=Y.*Y.*SCAL.*SCAL;
590Err=sum(temp);
591Err = max(sqrt(Err/N),1e-10);
592return
593
594
595function SDIRK_ErrorMsg(code,T,H,Ierr)
596Ierr = code;
597
598disp('\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~');
599disp('\nForced exit from Sdirk due to the following error:\n');
600
601switch (code)
602
603case -1
604        disp('--> Improper value for maximal no of steps');
605       
606case -2
607        disp('--> Selected SDIRK method not implemented');
608       
609case -3
610        printf('--> Hmin/Hmax/Hstart must be positive');
611       
612case -4
613        printf('--> FacMin/FacMax/FacRej must be positive');
614       
615case -5
616        printf('-> Improper tolerance values');
617case -6
618        disp('--> No of steps exceeds maximum bound');
619case -7
620        disp('--> Step size too small (T + H/10 = T) or H < Roundoff');
621case -8
622        printf('--> Matrix is repeatedly singular');
623
624end;
625return
626
627
628function [E1L,E1U,ISTATUS,ISING]=PrepareMatrix(N,H,T,Y,FJAC,SkipJac, SkipLU,Reject,ISTATUS,OdeJacobian1)
629global ZERO Njac  Ndec  Nsng
630global rkGamma
631
632
633ConsecutiveSng = 0;
634ISING =1;
635while (ISING ~= ZERO)
636    HGammaInv = 1/(H*rkGamma);
637    if(SkipJac == 0)
638        FJAC = OdeJacobian1(T,Y);
639        ISTATUS(Njac)=ISTATUS(Njac)+1;
640    end;
641    E=eye(N);
642    E=HGammaInv*E-FJAC;
643   
644    [E1L,E1U] = lu(E);
645    if(min(rank(E1L),rank(E1U))<N)
646        ISING = 1;
647    else
648        ISING = 0;
649    end;
650    ISTATUS(Ndec) =ISTATUS(Ndec)+1;
651    if(ISING ~= 0 )
652        ISTATUS(Nsng)=ISTATUS(Nsng)+1;
653        ConsecutiveSng=ConsecutiveSng+1;
654       
655        if(ConsecutiveSng >= 6)
656            return;
657        end
658        H = 0.5*H;
659        SkipJac =1;
660        SkipLU = 0;
661        Reject =1 ;
662    end;
663end
664return;
665
666
667function [RHS,ISTATUS] = SDIRK_Solve(H,N,E1L,E1U,RHS,ISTATUS)
668global  Nsol
669global  rkGamma
670
671HGammaInv = 1/(H*rkGamma);
672RHS = HGammaInv*RHS;
673size(RHS);
674size(E1L);
675size(E1U);
676temp = E1L\RHS;
677RHS = E1U\temp;
678ISTATUS(Nsol)=ISTATUS(Nsol)+1;
679return
680
681function Sdirk4a()
682%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*
683global ZERO ONE
684global  S4A  sdMethod rkS rkGamma rkA rkB rkELO rkBhat
685global rkC rkD rkE rkTheta rkAlpha
686sdMethod = S4A;
687
688%/* Number of stages */
689rkS = 5;
690
691%/* Method Coefficients */
692rkGamma = 0.2666666666666666666666666666666667;
693
694rkA(1,1) = 0.2666666666666666666666666666666667;
695rkA(1,2) = 0.5000000000000000000000000000000000;
696rkA(2,2) = 0.2666666666666666666666666666666667;
697rkA(1,3) = 0.3541539528432732316227461858529820;
698rkA(2,3) = -0.5415395284327323162274618585298197e-01;
699rkA(3,3) = 0.2666666666666666666666666666666667;
700rkA(1,4) = 0.8515494131138652076337791881433756e-01;
701rkA(2,4) = -0.6484332287891555171683963466229754e-01;
702rkA(3,4) = 0.7915325296404206392428857585141242e-01;
703rkA(4,4) = 0.2666666666666666666666666666666667;
704rkA(1,5) = 2.100115700566932777970612055999074;
705rkA(2,5) = -0.7677800284445976813343102185062276;
706rkA(3,5) = 2.399816361080026398094746205273880;
707rkA(4,5) = -2.998818699869028161397714709433394;
708rkA(5,5) = 0.2666666666666666666666666666666667;
709rkB(1) = 2.100115700566932777970612055999074;
710rkB(2) = -0.7677800284445976813343102185062276;
711rkB(3) = 2.399816361080026398094746205273880;
712rkB(4) = -2.998818699869028161397714709433394;
713rkB(5) = 0.2666666666666666666666666666666667;
714
715rkBhat(1) = 2.885264204387193942183851612883390;
716rkBhat(2) = -0.1458793482962771337341223443218041;
717rkBhat(3) = 2.390008682465139866479830743628554;
718rkBhat(4) = -4.129393538556056674929560012190140;
719rkBhat(5) = ZERO;
720
721rkC(1) = 0.2666666666666666666666666666666667;
722rkC(2) = 0.7666666666666666666666666666666667;
723rkC(3) = 0.5666666666666666666666666666666667;
724rkC(4) = 0.3661315380631796996374935266701191;
725rkC(5) = ONE;
726
727%/* Ynew = Yold + h*Sum_i {rkB_i*k_i} = Yold + Sum_i {rkD_i*Z_i} */
728rkD(1) = ZERO;
729rkD(2) = ZERO;
730rkD(3) = ZERO;
731rkD(4) = ZERO;
732rkD(5) = ONE;
733
734%/* Err = h * Sum_i {(rkB_i-rkBhat_i)*k_i} = Sum_i {rkE_i*Z_i} */
735rkE(1) = -0.6804000050475287124787034884002302;
736rkE(2) = 1.558961944525217193393931795738823;
737rkE(3) = -13.55893003128907927748632408763868;
738rkE(4) = 15.48522576958521253098585004571302;
739rkE(5) = ONE;
740
741%/* Local order of Err estimate */
742rkELO = 4;
743
744%/* h*Sum_j {rkA_ij*k_j} = Sum_j {rkTheta_ij*Z_j} */
745rkTheta(1,2) = 1.875000000000000000000000000000000;
746rkTheta(1,3) = 1.708847304091539528432732316227462;
747rkTheta(2,3) = -0.2030773231622746185852981969486824;
748rkTheta(1,4) = 0.2680325578937783958847157206823118;
749rkTheta(2,4) = -0.1828840955527181631794050728644549;
750rkTheta(3,4) = 0.2968246986151577397160821594427966;
751rkTheta(1,5) = 0.9096171815241460655379433581446771;
752rkTheta(2,5) = -3.108254967778352416114774430509465;
753rkTheta(3,5) = 12.33727431701306195581826123274001;
754rkTheta(4,5) = -11.24557012450885560524143016037523;
755
756%/* Starting value for Newton iterations: Z_i^0 = Sum_j {rkAlpha_ij*Z_j} */
757rkAlpha(1,2) = 2.875000000000000000000000000000000;
758rkAlpha(1,3) = 0.8500000000000000000000000000000000;
759rkAlpha(2,3) = 0.4434782608695652173913043478260870;
760rkAlpha(1,4) = 0.7352046091658870564637910527807370;
761rkAlpha(2,4) = -0.9525565003057343527941920657462074e-01;
762rkAlpha(3,4) = 0.4290111305453813852259481840631738;
763rkAlpha(1,5) = -16.10898993405067684831655675112808;
764rkAlpha(2,5) = 6.559571569643355712998131800797873;
765rkAlpha(3,5) = -15.90772144271326504260996815012482;
766rkAlpha(4,5) = 25.34908987169226073668861694892683;
767
768rkELO = 4.0;
769
770% /* end Sdirk4a */
771
772%/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
773function Sdirk4b()
774%/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
775global ZERO ONE
776global  S4B sdMethod rkS rkGamma rkA rkB rkELO rkBhat
777global rkC rkD rkE rkTheta rkAlpha
778
779sdMethod = S4B;
780
781%/* Number of stages */
782rkS = 5;
783
784%/* Method coefficients */
785rkGamma = 0.25;
786
787rkA(1,1) = 0.25;
788rkA(1,2) = 0.5;
789rkA(2,2) = 0.25;
790rkA(1,3) = 0.34;
791rkA(2,3) = -0.40e-01;
792rkA(3,3) = 0.25;
793rkA(1,4) = 0.2727941176470588235294117647058824;
794rkA(2,4) = -0.5036764705882352941176470588235294e-01;
795rkA(3,4) = 0.2757352941176470588235294117647059e-01;
796rkA(4,4) = 0.25;
797rkA(1,5) = 1.041666666666666666666666666666667;
798rkA(2,5) = -1.020833333333333333333333333333333;
799rkA(3,5) = 7.812500000000000000000000000000000;
800rkA(4,5) = -7.083333333333333333333333333333333;
801rkA(5,5) = 0.25;
802
803rkB(1) = 1.041666666666666666666666666666667;
804rkB(2) = -1.020833333333333333333333333333333;
805rkB(3) = 7.812500000000000000000000000000000;
806rkB(4) = -7.083333333333333333333333333333333;
807rkB(5) = 0.250000000000000000000000000000000;
808
809rkBhat(1) = 1.069791666666666666666666666666667;
810rkBhat(2) = -0.894270833333333333333333333333333;
811rkBhat(3) = 7.695312500000000000000000000000000;
812rkBhat(4) = -7.083333333333333333333333333333333;
813rkBhat(5) = 0.212500000000000000000000000000000;
814
815rkC(1) = 0.25;
816rkC(2) = 0.75;
817rkC(3) = 0.55;
818rkC(4) = 0.5;
819rkC(5) = ONE;
820
821%/* Ynew = Yold + h*Sum_i {rkB_i*k_i} = Yold + Sum_i {rkD_i*Z_i} */
822rkD(1) = ZERO;
823rkD(2) = ZERO;
824rkD(3) = ZERO;
825rkD(4) = ZERO;
826rkD(5) = ONE;
827
828%/* Err = h * Sum_i {(rkB_i-rkBhat_i)*k_i} = Sum_i {rkE_i*Z_i} */
829rkE(1) = 0.5750;
830rkE(2) = 0.2125;
831rkE(3) = -4.6875;
832rkE(4) = 4.2500;
833rkE(5) = 0.1500;
834
835%/* Local order of Err estimate */
836rkELO = 4;
837
838%/* h*Sum_j {rkA_ij*k_j} = Sum_j {rkTheta_ij*Z_j} */
839rkTheta(1,2) = 2.0;
840rkTheta(1,3) = 1.680000000000000000000000000000000;
841rkTheta(2,3) = -0.1600000000000000000000000000000000;
842rkTheta(1,4) = 1.308823529411764705882352941176471;
843rkTheta(2,4) = -0.1838235294117647058823529411764706;
844rkTheta(3,4) = 0.1102941176470588235294117647058824;
845rkTheta(1,5) = -3.083333333333333333333333333333333;
846rkTheta(2,5) = -4.291666666666666666666666666666667;
847rkTheta(3,5) = 34.37500000000000000000000000000000;
848rkTheta(4,5) = -28.3333333333333333333333333333;
849
850%/* Starting value for Newton iterations: Z_i^0 = Sum_j {rkAlpha_ij*Z_j} */
851rkAlpha(1,2) = 3.0;
852rkAlpha(1,3) = 0.8800000000000000000000000000000000;
853rkAlpha(2,3) = 0.4400000000000000000000000000000000;
854rkAlpha(1,4) = 0.1666666666666666666666666666666667;
855rkAlpha(2,4) = -0.8333333333333333333333333333333333e-01;
856rkAlpha(3,4) = 0.9469696969696969696969696969696970;
857rkAlpha(1,5) = -6.0;
858rkAlpha(2,5) = 9.0;
859rkAlpha(3,5) = -56.81818181818181818181818181818182;
860rkAlpha(4,5) = 54.0;
861
862 %/* end Sdirk4b */
863
864%/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
865function Sdirk2a()
866%/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
867global ZERO ONE
868global  S2A  sdMethod rkS rkGamma rkA rkB rkELO rkBhat
869global rkC rkD rkE rkTheta rkAlpha
870
871sdMethod = S2A;
872
873%/* ~~~> Number of stages    */
874rkS = 2;
875
876%/* ~~~> Method coefficients    */
877rkGamma = 0.2928932188134524755991556378951510;
878rkA(1,1) = 0.2928932188134524755991556378951510;
879rkA(1,2) = 0.7071067811865475244008443621048490;
880rkA(2,2) = 0.2928932188134524755991556378951510;
881rkB(1) = 0.7071067811865475244008443621048490;
882rkB(2) = 0.2928932188134524755991556378951510;
883rkBhat(1) = 0.6666666666666666666666666666666667;
884rkBhat(2) = 0.3333333333333333333333333333333333;
885rkC(1) = 0.292893218813452475599155637895151;
886rkC(2) = ONE;
887
888%/* ~~~> Ynew = Yold + h*Sum_i {rkB_i*k_i} = Yold + Sum_i {rkD_i*Z_i}    */
889rkD(1) = ZERO;
890rkD(2) = ONE;
891
892%/* ~~~> Err = h * Sum_i {(rkB_i-rkBhat_i)*k_i} = Sum_i {rkE_i*Z_i}    */
893rkE(1) = 0.4714045207910316829338962414032326;
894rkE(2) = -0.1380711874576983496005629080698993;
895
896%/* ~~~> Local order of Err estimate    */
897rkELO = 2;
898
899%/* ~~~> h*Sum_j {rkA_ij*k_j} = Sum_j {rkTheta_ij*Z_j}    */
900rkTheta(1,2) = 2.414213562373095048801688724209698;
901
902%/* ~~~> Starting value for Newton iterations */
903rkAlpha(1,2) = 3.414213562373095048801688724209698;
904
905% /* end Sdirk2a */
906
907%/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
908function Sdirk2b()
909%/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
910global ZERO ONE
911global  S2B  sdMethod rkS rkGamma rkA rkB rkELO rkBhat
912global rkC rkD rkE rkTheta rkAlpha
913
914sdMethod = S2B;
915
916%/* ~~~> Number of stages    */
917rkS      = 2;
918
919%/* ~~~> Method coefficients    */
920rkGamma = 1.707106781186547524400844362104849;
921rkA(1,1) = 1.707106781186547524400844362104849;
922rkA(1,2) = -0.707106781186547524400844362104849;
923rkA(2,2) = 1.707106781186547524400844362104849;
924rkB(1) = -0.707106781186547524400844362104849;
925rkB(2) = 1.707106781186547524400844362104849;
926rkBhat(1) = 0.6666666666666666666666666666666667;
927rkBhat(2) = 0.3333333333333333333333333333333333;
928rkC(1) = 1.707106781186547524400844362104849;
929rkC(2) = ONE;
930
931%/* ~~~> Ynew = Yold + h*Sum_i {rkB_i*k_i} = Yold + Sum_i {rkD_i*Z_i}    */
932rkD(1) = ZERO;
933rkD(2) = ONE;
934
935%/* ~~~> Err = h * Sum_i {(rkB_i-rkBhat_i)*k_i} = Sum_i {rkE_i*Z_i}    */
936rkE(1) = -0.4714045207910316829338962414032326;
937rkE(2) = 0.8047378541243650162672295747365659;
938
939%/* ~~~> Local order of Err estimate    */
940rkELO = 2;
941
942%/* ~~~> h*Sum_j {rkA_ij*k_j} = Sum_j {rkTheta_ij*Z_j}    */
943rkTheta(1,2) = -0.414213562373095048801688724209698;
944
945%/* ~~~> Starting value for Newton iterations */
946rkAlpha(1,2) = 0.5857864376269049511983112757903019;
947
948% /* end Sdirk2b */
949
950%/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
951function Sdirk3a()
952%/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
953global ZERO ONE
954global  S3A  sdMethod rkS rkGamma rkA rkB rkELO rkBhat
955global rkC rkD rkE rkTheta rkAlpha
956
957sdMethod = S3A;
958
959%/* ~~~> Number of stages    */
960rkS = 3;
961
962%/* ~~~> Method coefficients    */
963rkGamma = 0.2113248654051871177454256097490213;
964rkA(1,1) = 0.2113248654051871177454256097490213;
965rkA(1,2) = 0.2113248654051871177454256097490213;
966rkA(2,2) = 0.2113248654051871177454256097490213;
967rkA(1,3) = 0.2113248654051871177454256097490213;
968rkA(2,3) = 0.5773502691896257645091487805019573;
969rkA(3,3) = 0.2113248654051871177454256097490213;
970rkB(1) = 0.2113248654051871177454256097490213;
971rkB(2) = 0.5773502691896257645091487805019573;
972rkB(3) = 0.2113248654051871177454256097490213;
973rkBhat(1)= 0.2113248654051871177454256097490213;
974rkBhat(2)= 0.6477918909913548037576239837516312;
975rkBhat(3)= 0.1408832436034580784969504064993475;
976rkC(1) = 0.2113248654051871177454256097490213;
977rkC(2) = 0.4226497308103742354908512194980427;
978rkC(3) = ONE;
979
980%/* ~~~> Ynew = Yold + h*Sum_i {rkB_i*k_i} = Yold + Sum_i {rkD_i*Z_i}    */
981rkD(1) = ZERO;
982rkD(2) = ZERO;
983rkD(3) = ONE;
984
985%/* ~~~> Err = h * Sum_i {(rkB_i-rkBhat_i)*k_i} = Sum_i {rkE_i*Z_i}    */
986rkE(1) = 0.9106836025229590978424821138352906;
987rkE(2) = -1.244016935856292431175815447168624;
988rkE(3) = 0.3333333333333333333333333333333333;
989
990%/* ~~~> Local order of Err estimate    */
991rkELO    = 2;
992
993%/* ~~~> h*Sum_j {rkA_ij*k_j} = Sum_j {rkTheta_ij*Z_j}    */
994rkTheta(1,2) =  ONE;
995rkTheta(1,3) = -1.732050807568877293527446341505872;
996rkTheta(2,3) = 2.732050807568877293527446341505872;
997
998%/* ~~~> Starting value for Newton iterations */
999rkAlpha(1,2) = 2.0;
1000rkAlpha(1,3) = -12.92820323027550917410978536602349;
1001rkAlpha(2,3) = 8.83012701892219323381861585376468;
1002
1003% /* end Sdirk3a */   
1004
1005
1006     
Note: See TracBrowser for help on using the repository browser.