[2696] | 1 | function [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 | %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 109 | global Nfun Njac Nstp Nacc Nrej Ndec Nsol Nsng Ntexit Nhexit Nhnew |
---|
| 110 | % Parse ODE options |
---|
| 111 | AbsTol = odeget(Options, 'AbsTol'); |
---|
| 112 | if isempty(AbsTol) |
---|
| 113 | AbsTol = 1.0e-3; |
---|
| 114 | end |
---|
| 115 | Hstart = odeget(Options, 'InitialStep'); |
---|
| 116 | if isempty(Hstart) |
---|
| 117 | Hstart = 0; |
---|
| 118 | end |
---|
| 119 | Jacobian = odeget(Options, 'Jacobian'); |
---|
| 120 | if isempty(Jacobian) |
---|
| 121 | error('A Jacobian function is required.'); |
---|
| 122 | end |
---|
| 123 | Hmax = odeget(Options, 'MaxStep'); |
---|
| 124 | if isempty(Hmax) |
---|
| 125 | Hmax = 0; |
---|
| 126 | end |
---|
| 127 | RelTol = odeget(Options, 'RelTol'); |
---|
| 128 | if isempty(RelTol) |
---|
| 129 | RelTol = 1.0e-4; |
---|
| 130 | end |
---|
| 131 | % Initialize statistics |
---|
| 132 | Nfun=1; Njac=2; Nstp=3; Nacc=4; Nrej=5; Ndec=6; Nsol=7; Nsng=8; |
---|
| 133 | Ntexit=1; Nhexit=2; Nhnew=3; |
---|
| 134 | RSTATUS = zeros(20,1); |
---|
| 135 | ISTATUS = zeros(20,1); |
---|
| 136 | |
---|
| 137 | % Get problem size |
---|
| 138 | steps = length(Tspan); |
---|
| 139 | N = max(size(Y0)); |
---|
| 140 | |
---|
| 141 | % Initialize tolerances |
---|
| 142 | ATOL = ones(N,1)*AbsTol; |
---|
| 143 | RTOL = ones(N,1)*RelTol; |
---|
| 144 | |
---|
| 145 | % Initialize step |
---|
| 146 | RCNTRL(2) = max(0, Hmax); |
---|
| 147 | RCNTRL(3) = max(0, Hstart); |
---|
| 148 | |
---|
| 149 | % Integrate |
---|
| 150 | Y = zeros(N,steps); |
---|
| 151 | T = zeros(steps,1); |
---|
| 152 | Y(:,1) = Y0; |
---|
| 153 | T(1) = Tspan(1); |
---|
| 154 | t=cputime; |
---|
| 155 | for 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 | |
---|
| 164 | end |
---|
| 165 | cputime-t |
---|
| 166 | |
---|
| 167 | |
---|
| 168 | |
---|
| 169 | function [T,Y,Ierr,ISTATUS,RSTATUS] = Sdirk_Integrate(N,Y0,TIN, TOUT,OdeFunction1,... |
---|
| 170 | OdeJacobian1,ATOL,RTOL,RCNTRL,ICNTRL,ISTATUS,RSTATUS) |
---|
| 171 | global ZERO ONE Nfun Njac Nstp Nacc Nrej Ndec Nsol Nsng Ntexit Nhexit Nhnew |
---|
| 172 | global Smax S2A S2B S3A S4A S4B |
---|
| 173 | global rkTheta rkAlpha |
---|
| 174 | ZERO = 0; ONE=1.0; Nfun = 1; Njac = 2; Nstp = 3; Nacc = 4; Nrej = 5;Ndec=6; |
---|
| 175 | Nsol = 7; Nsng = 8; Ntexit = 1; Nhexit = 2; Nhnew = 3; Smax = 5; S2A = 1; |
---|
| 176 | S2B = 2; S3A = 3; S4A = 4; S4B = 5; |
---|
| 177 | rkTheta = zeros(Smax); |
---|
| 178 | rkAlpha = zeros(Smax); |
---|
| 179 | |
---|
| 180 | %Local Variables |
---|
| 181 | ISTATUS = zeros(20,1); |
---|
| 182 | RSTATUS = zeros(20,1); |
---|
| 183 | |
---|
| 184 | % Fine-Tune the integrator |
---|
| 185 | T1 = TIN; |
---|
| 186 | T2 = TOUT; |
---|
| 187 | RTOL=ones(size(Y0)); |
---|
| 188 | RTOL=1e-3*RTOL; |
---|
| 189 | ATOL=RTOL; |
---|
| 190 | |
---|
| 191 | [T,Y,Ierr,RSTATUS,ISTATUS] = SDIRK(N, T1,T2,Y0,RTOL,ATOL,RCNTRL,ICNTRL,RSTATUS,ISTATUS,OdeFunction1,OdeJacobian1); |
---|
| 192 | |
---|
| 193 | if Ierr < 0 |
---|
| 194 | disp(['SDIRK:Unsuccessful exit at ',num2str(TIN), ' ',num2str(Ierr)]); |
---|
| 195 | end |
---|
| 196 | return; |
---|
| 197 | |
---|
| 198 | |
---|
| 199 | |
---|
| 200 | |
---|
| 201 | function [T, Y, Ierr,RSTATUS,ISTATUS] = SDIRK(N,T1,T2,Y,AbsTol,RelTol,... |
---|
| 202 | RCNTRL,ICNTRL,RSTATUS,ISTATUS,OdeFunction1,OdeJacobian1) |
---|
| 203 | |
---|
| 204 | global ZERO ONE |
---|
| 205 | |
---|
| 206 | |
---|
| 207 | Max_no_steps = 0;Hmin = 0;Hmax = 0; Hstart = 0; Roundoff = eps; FacMin = 0; |
---|
| 208 | FacMax = 0; FacSafe = 0; FacRej = 0; |
---|
| 209 | if(ICNTRL(1)==0) |
---|
| 210 | ITOL = 1; |
---|
| 211 | else |
---|
| 212 | ITOL = 0; |
---|
| 213 | end |
---|
| 214 | Ierr = 0; |
---|
| 215 | %~~~> Method Selection |
---|
| 216 | switch (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(); |
---|
| 232 | end |
---|
| 233 | |
---|
| 234 | if (ICNTRL(3)==ZERO) |
---|
| 235 | Max_no_steps = 200000; |
---|
| 236 | end |
---|
| 237 | if(ICNTRL(3) > ZERO) |
---|
| 238 | Max_no_steps = ICNTRL(3); |
---|
| 239 | end |
---|
| 240 | |
---|
| 241 | if(ICNTRL(3) < ZERO) |
---|
| 242 | disp(['User-Selected ICNTRL(3)=',num2str(ICNTRL(3))]); |
---|
| 243 | SDIRK_ErrorMsg(-1,T1,ZERO,Ierr); |
---|
| 244 | end; |
---|
| 245 | %~~>StartNewton Extrapolate for starting values of Newton iterations |
---|
| 246 | if(ICNTRL(5)==ZERO) |
---|
| 247 | StartNewton = 1; |
---|
| 248 | else |
---|
| 249 | StartNewton = 0; |
---|
| 250 | end; |
---|
| 251 | |
---|
| 252 | if(ICNTRL(4)>=ZERO) |
---|
| 253 | if(ICNTRL(4)==ZERO) |
---|
| 254 | NewtonMaxit = 8; |
---|
| 255 | else |
---|
| 256 | NewtonMaxit = ICNTRL(4); |
---|
| 257 | end |
---|
| 258 | else |
---|
| 259 | disp(['User-Selected NewtonMaxit ICNTRL(4)=',num2str(ICNTRL(4))]); |
---|
| 260 | SDIRK_ErrorMsg(-2,T1,ZERO,Ierr); |
---|
| 261 | end; |
---|
| 262 | |
---|
| 263 | %~~> Lower bound on the step size (Positive value) |
---|
| 264 | |
---|
| 265 | if(RCNTRL(1)==ZERO) |
---|
| 266 | Hmin = ZERO; |
---|
| 267 | end |
---|
| 268 | if (RCNTRL(1)>ZERO) |
---|
| 269 | Hmin = RCNTRL(1); |
---|
| 270 | end |
---|
| 271 | if(RCNTRL(1)<ZERO) |
---|
| 272 | disp(['User-Selected Hmin, RCNTRL(1)=',num2str(RCNTRL(1))]); |
---|
| 273 | SDIRK_ErrorMsg(-3,T1,ZERO,Ierr); |
---|
| 274 | end; |
---|
| 275 | % Upper bound on the step size : (Positive value) |
---|
| 276 | |
---|
| 277 | if(RCNTRL(2)==ZERO) |
---|
| 278 | Hmax = abs(T2-T1); |
---|
| 279 | end |
---|
| 280 | if(RCNTRL(2)>ZERO) |
---|
| 281 | Hmax = min(abs(T2-T1),RCNTRL(2)); |
---|
| 282 | end; |
---|
| 283 | if(RCNTRL(2)<ZERO) |
---|
| 284 | disp(['User-Selected Hmax, RCNTRL(2)=',num2str(RCNTRL(2))]); |
---|
| 285 | SDIRK_ErrorMsg(-3,T1,ZERO,Ierr); |
---|
| 286 | end; |
---|
| 287 | %~~> Starting step size |
---|
| 288 | |
---|
| 289 | if(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 |
---|
| 295 | else |
---|
| 296 | disp(['User-Selected Hstart, RCNTRL(3)=',num2str(RCNTRL(3))]); |
---|
| 297 | SDIRK_ErrorMsg(-3,T1,ZERO,Ierr); |
---|
| 298 | end; |
---|
| 299 | %~~> Stepsize can be changed such that FacMin < Hnew/Hexit < FacMax |
---|
| 300 | if(RCNTRL(4) >= ZERO) |
---|
| 301 | if(RCNTRL(4) == ZERO) |
---|
| 302 | FacMin = 0.2; |
---|
| 303 | else |
---|
| 304 | FacMin = RCNTRL(4); |
---|
| 305 | end |
---|
| 306 | else |
---|
| 307 | disp(['User-Selected FacMin, RCNTRL(4)=',num2str(RCNTRL(4))]); |
---|
| 308 | SDIRK_ErrorMsg(-4,T1,ZERO,Ierr); |
---|
| 309 | end |
---|
| 310 | |
---|
| 311 | |
---|
| 312 | if (RCNTRL(5)>=ZERO) |
---|
| 313 | if(RCNTRL(5)==ZERO) |
---|
| 314 | FacMax = 10.0; |
---|
| 315 | else |
---|
| 316 | FacMax = RCNTRL(5); |
---|
| 317 | end |
---|
| 318 | else |
---|
| 319 | disp(['User-Selected FacMax, RCNTRL(5)=',num2str(RCNTRL(5))]); |
---|
| 320 | SDIRK_ErrorMsg(-4,T1,ZERO,Ierr); |
---|
| 321 | end |
---|
| 322 | |
---|
| 323 | %~~> FacRej Factor to decrease step after 2 successive rejections |
---|
| 324 | |
---|
| 325 | if(RCNTRL(6)>=ZERO) |
---|
| 326 | if(RCNTRL(6) == ZERO) |
---|
| 327 | FacRej = 0.1; |
---|
| 328 | else |
---|
| 329 | FacRej = RCNTRL(6); |
---|
| 330 | end |
---|
| 331 | else |
---|
| 332 | disp(['User-Selected FacRej, RCNTRL(6)=',num2str(RCNTRL(6))]); |
---|
| 333 | SDIRK_ErrorMsg(-4,T1,ZERO,Ierr); |
---|
| 334 | end; |
---|
| 335 | |
---|
| 336 | %~~> FacSafe Safety Factor in the computation of new step size |
---|
| 337 | |
---|
| 338 | if(RCNTRL(7)>=ZERO) |
---|
| 339 | if(RCNTRL(7)==ZERO) |
---|
| 340 | FacSafe = 0.9; |
---|
| 341 | else |
---|
| 342 | FacSafe = RCNTRL(7); |
---|
| 343 | end |
---|
| 344 | else |
---|
| 345 | disp(['User-Selected FacSafe, RCNTRL(7)=',num2str(RCNTRL(7))]); |
---|
| 346 | SDIRK_ErrorMsg(-4,T1,ZERO,Ierr); |
---|
| 347 | end |
---|
| 348 | |
---|
| 349 | %~~> ThetaMin decides whether the Jacobian should be recomputed |
---|
| 350 | |
---|
| 351 | if(RCNTRL(8) == ZERO) |
---|
| 352 | ThetaMin = 1.0E-3; |
---|
| 353 | else |
---|
| 354 | ThetaMin = RCNTRL(8); |
---|
| 355 | end; |
---|
| 356 | |
---|
| 357 | %~~> Stopping Criteria for Newtons method |
---|
| 358 | |
---|
| 359 | if(RCNTRL(9) == ZERO) |
---|
| 360 | NewtonTol = 3.0e-4; |
---|
| 361 | else |
---|
| 362 | NewtonTol = RCNTRL(9); |
---|
| 363 | end; |
---|
| 364 | |
---|
| 365 | %~~> Qmin, Qmax : IF Qmin <Hnew/Hold<Qmax STEP SIZE = constant |
---|
| 366 | |
---|
| 367 | if(RCNTRL(10) == ZERO) |
---|
| 368 | Qmin = ONE; |
---|
| 369 | else |
---|
| 370 | Qmin = RCNTRL(10); |
---|
| 371 | end |
---|
| 372 | |
---|
| 373 | if(RCNTRL(11) == ZERO) |
---|
| 374 | Qmax = 1.2; |
---|
| 375 | else |
---|
| 376 | Qmax = RCNTRL(11); |
---|
| 377 | end |
---|
| 378 | |
---|
| 379 | %~~> Check if tolerances are reasonable |
---|
| 380 | if(ITOL == ZERO) |
---|
| 381 | if((AbsTol(1)<=ZERO) || RelTol(1) <=10*Roundoff) |
---|
| 382 | SDIRK_ErrorMsg(-5,T1,ZERO,Ierr); |
---|
| 383 | end; |
---|
| 384 | else |
---|
| 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 |
---|
| 390 | end |
---|
| 391 | if(Ierr < 0) |
---|
| 392 | return; |
---|
| 393 | end |
---|
| 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 | |
---|
| 404 | return; |
---|
| 405 | |
---|
| 406 | function[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) |
---|
| 410 | global ZERO ONE Nfun Nstp Nacc Nrej Ntexit Nhexit Nhnew |
---|
| 411 | global rkS rkGamma rkELO |
---|
| 412 | global rkC rkD rkE rkTheta rkAlpha |
---|
| 413 | T =T1; |
---|
| 414 | Tdirection = sign(T2-T1); |
---|
| 415 | H = max(abs(Hmin),abs(Hstart)); |
---|
| 416 | if(H<=10.0*Roundoff) |
---|
| 417 | H = 1.0e-6; |
---|
| 418 | end; |
---|
| 419 | H = min(H,Hmax); |
---|
| 420 | H=sign(Tdirection)*H; |
---|
| 421 | SkipLU = 0; |
---|
| 422 | SkipJac = 0; |
---|
| 423 | Reject = 0; |
---|
| 424 | FirstStep = 1; |
---|
| 425 | CycleTloop = 0; |
---|
| 426 | SCAL = SDIRK_ErrorScale(N,ITOL,AbsTol,RelTol,Y); |
---|
| 427 | FJAC = zeros(N); |
---|
| 428 | while((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 |
---|
| 567 | end |
---|
| 568 | |
---|
| 569 | Ierr =1; |
---|
| 570 | return; |
---|
| 571 | |
---|
| 572 | |
---|
| 573 | |
---|
| 574 | |
---|
| 575 | function [SCAL] = SDIRK_ErrorScale(N,ITOL,AbsTol,RelTol,Y) |
---|
| 576 | ZERO =0; |
---|
| 577 | if(ITOL==ZERO) |
---|
| 578 | temp = RelTol(1)*Y; |
---|
| 579 | temp =(AbsTol(1)+abs(temp)); |
---|
| 580 | SCAL=1./temp; |
---|
| 581 | else |
---|
| 582 | SCAL=1./(AbsTol+RelTol.*abs(Y)); |
---|
| 583 | end |
---|
| 584 | return; |
---|
| 585 | |
---|
| 586 | |
---|
| 587 | function [Err] = SDIRK_ErrorNorm(N,Y,SCAL) |
---|
| 588 | ZERO=0; |
---|
| 589 | temp=Y.*Y.*SCAL.*SCAL; |
---|
| 590 | Err=sum(temp); |
---|
| 591 | Err = max(sqrt(Err/N),1e-10); |
---|
| 592 | return |
---|
| 593 | |
---|
| 594 | |
---|
| 595 | function SDIRK_ErrorMsg(code,T,H,Ierr) |
---|
| 596 | Ierr = code; |
---|
| 597 | |
---|
| 598 | disp('\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'); |
---|
| 599 | disp('\nForced exit from Sdirk due to the following error:\n'); |
---|
| 600 | |
---|
| 601 | switch (code) |
---|
| 602 | |
---|
| 603 | case -1 |
---|
| 604 | disp('--> Improper value for maximal no of steps'); |
---|
| 605 | |
---|
| 606 | case -2 |
---|
| 607 | disp('--> Selected SDIRK method not implemented'); |
---|
| 608 | |
---|
| 609 | case -3 |
---|
| 610 | printf('--> Hmin/Hmax/Hstart must be positive'); |
---|
| 611 | |
---|
| 612 | case -4 |
---|
| 613 | printf('--> FacMin/FacMax/FacRej must be positive'); |
---|
| 614 | |
---|
| 615 | case -5 |
---|
| 616 | printf('-> Improper tolerance values'); |
---|
| 617 | case -6 |
---|
| 618 | disp('--> No of steps exceeds maximum bound'); |
---|
| 619 | case -7 |
---|
| 620 | disp('--> Step size too small (T + H/10 = T) or H < Roundoff'); |
---|
| 621 | case -8 |
---|
| 622 | printf('--> Matrix is repeatedly singular'); |
---|
| 623 | |
---|
| 624 | end; |
---|
| 625 | return |
---|
| 626 | |
---|
| 627 | |
---|
| 628 | function [E1L,E1U,ISTATUS,ISING]=PrepareMatrix(N,H,T,Y,FJAC,SkipJac, SkipLU,Reject,ISTATUS,OdeJacobian1) |
---|
| 629 | global ZERO Njac Ndec Nsng |
---|
| 630 | global rkGamma |
---|
| 631 | |
---|
| 632 | |
---|
| 633 | ConsecutiveSng = 0; |
---|
| 634 | ISING =1; |
---|
| 635 | while (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; |
---|
| 663 | end |
---|
| 664 | return; |
---|
| 665 | |
---|
| 666 | |
---|
| 667 | function [RHS,ISTATUS] = SDIRK_Solve(H,N,E1L,E1U,RHS,ISTATUS) |
---|
| 668 | global Nsol |
---|
| 669 | global rkGamma |
---|
| 670 | |
---|
| 671 | HGammaInv = 1/(H*rkGamma); |
---|
| 672 | RHS = HGammaInv*RHS; |
---|
| 673 | size(RHS); |
---|
| 674 | size(E1L); |
---|
| 675 | size(E1U); |
---|
| 676 | temp = E1L\RHS; |
---|
| 677 | RHS = E1U\temp; |
---|
| 678 | ISTATUS(Nsol)=ISTATUS(Nsol)+1; |
---|
| 679 | return |
---|
| 680 | |
---|
| 681 | function Sdirk4a() |
---|
| 682 | %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~* |
---|
| 683 | global ZERO ONE |
---|
| 684 | global S4A sdMethod rkS rkGamma rkA rkB rkELO rkBhat |
---|
| 685 | global rkC rkD rkE rkTheta rkAlpha |
---|
| 686 | sdMethod = S4A; |
---|
| 687 | |
---|
| 688 | %/* Number of stages */ |
---|
| 689 | rkS = 5; |
---|
| 690 | |
---|
| 691 | %/* Method Coefficients */ |
---|
| 692 | rkGamma = 0.2666666666666666666666666666666667; |
---|
| 693 | |
---|
| 694 | rkA(1,1) = 0.2666666666666666666666666666666667; |
---|
| 695 | rkA(1,2) = 0.5000000000000000000000000000000000; |
---|
| 696 | rkA(2,2) = 0.2666666666666666666666666666666667; |
---|
| 697 | rkA(1,3) = 0.3541539528432732316227461858529820; |
---|
| 698 | rkA(2,3) = -0.5415395284327323162274618585298197e-01; |
---|
| 699 | rkA(3,3) = 0.2666666666666666666666666666666667; |
---|
| 700 | rkA(1,4) = 0.8515494131138652076337791881433756e-01; |
---|
| 701 | rkA(2,4) = -0.6484332287891555171683963466229754e-01; |
---|
| 702 | rkA(3,4) = 0.7915325296404206392428857585141242e-01; |
---|
| 703 | rkA(4,4) = 0.2666666666666666666666666666666667; |
---|
| 704 | rkA(1,5) = 2.100115700566932777970612055999074; |
---|
| 705 | rkA(2,5) = -0.7677800284445976813343102185062276; |
---|
| 706 | rkA(3,5) = 2.399816361080026398094746205273880; |
---|
| 707 | rkA(4,5) = -2.998818699869028161397714709433394; |
---|
| 708 | rkA(5,5) = 0.2666666666666666666666666666666667; |
---|
| 709 | rkB(1) = 2.100115700566932777970612055999074; |
---|
| 710 | rkB(2) = -0.7677800284445976813343102185062276; |
---|
| 711 | rkB(3) = 2.399816361080026398094746205273880; |
---|
| 712 | rkB(4) = -2.998818699869028161397714709433394; |
---|
| 713 | rkB(5) = 0.2666666666666666666666666666666667; |
---|
| 714 | |
---|
| 715 | rkBhat(1) = 2.885264204387193942183851612883390; |
---|
| 716 | rkBhat(2) = -0.1458793482962771337341223443218041; |
---|
| 717 | rkBhat(3) = 2.390008682465139866479830743628554; |
---|
| 718 | rkBhat(4) = -4.129393538556056674929560012190140; |
---|
| 719 | rkBhat(5) = ZERO; |
---|
| 720 | |
---|
| 721 | rkC(1) = 0.2666666666666666666666666666666667; |
---|
| 722 | rkC(2) = 0.7666666666666666666666666666666667; |
---|
| 723 | rkC(3) = 0.5666666666666666666666666666666667; |
---|
| 724 | rkC(4) = 0.3661315380631796996374935266701191; |
---|
| 725 | rkC(5) = ONE; |
---|
| 726 | |
---|
| 727 | %/* Ynew = Yold + h*Sum_i {rkB_i*k_i} = Yold + Sum_i {rkD_i*Z_i} */ |
---|
| 728 | rkD(1) = ZERO; |
---|
| 729 | rkD(2) = ZERO; |
---|
| 730 | rkD(3) = ZERO; |
---|
| 731 | rkD(4) = ZERO; |
---|
| 732 | rkD(5) = ONE; |
---|
| 733 | |
---|
| 734 | %/* Err = h * Sum_i {(rkB_i-rkBhat_i)*k_i} = Sum_i {rkE_i*Z_i} */ |
---|
| 735 | rkE(1) = -0.6804000050475287124787034884002302; |
---|
| 736 | rkE(2) = 1.558961944525217193393931795738823; |
---|
| 737 | rkE(3) = -13.55893003128907927748632408763868; |
---|
| 738 | rkE(4) = 15.48522576958521253098585004571302; |
---|
| 739 | rkE(5) = ONE; |
---|
| 740 | |
---|
| 741 | %/* Local order of Err estimate */ |
---|
| 742 | rkELO = 4; |
---|
| 743 | |
---|
| 744 | %/* h*Sum_j {rkA_ij*k_j} = Sum_j {rkTheta_ij*Z_j} */ |
---|
| 745 | rkTheta(1,2) = 1.875000000000000000000000000000000; |
---|
| 746 | rkTheta(1,3) = 1.708847304091539528432732316227462; |
---|
| 747 | rkTheta(2,3) = -0.2030773231622746185852981969486824; |
---|
| 748 | rkTheta(1,4) = 0.2680325578937783958847157206823118; |
---|
| 749 | rkTheta(2,4) = -0.1828840955527181631794050728644549; |
---|
| 750 | rkTheta(3,4) = 0.2968246986151577397160821594427966; |
---|
| 751 | rkTheta(1,5) = 0.9096171815241460655379433581446771; |
---|
| 752 | rkTheta(2,5) = -3.108254967778352416114774430509465; |
---|
| 753 | rkTheta(3,5) = 12.33727431701306195581826123274001; |
---|
| 754 | rkTheta(4,5) = -11.24557012450885560524143016037523; |
---|
| 755 | |
---|
| 756 | %/* Starting value for Newton iterations: Z_i^0 = Sum_j {rkAlpha_ij*Z_j} */ |
---|
| 757 | rkAlpha(1,2) = 2.875000000000000000000000000000000; |
---|
| 758 | rkAlpha(1,3) = 0.8500000000000000000000000000000000; |
---|
| 759 | rkAlpha(2,3) = 0.4434782608695652173913043478260870; |
---|
| 760 | rkAlpha(1,4) = 0.7352046091658870564637910527807370; |
---|
| 761 | rkAlpha(2,4) = -0.9525565003057343527941920657462074e-01; |
---|
| 762 | rkAlpha(3,4) = 0.4290111305453813852259481840631738; |
---|
| 763 | rkAlpha(1,5) = -16.10898993405067684831655675112808; |
---|
| 764 | rkAlpha(2,5) = 6.559571569643355712998131800797873; |
---|
| 765 | rkAlpha(3,5) = -15.90772144271326504260996815012482; |
---|
| 766 | rkAlpha(4,5) = 25.34908987169226073668861694892683; |
---|
| 767 | |
---|
| 768 | rkELO = 4.0; |
---|
| 769 | |
---|
| 770 | % /* end Sdirk4a */ |
---|
| 771 | |
---|
| 772 | %/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 773 | function Sdirk4b() |
---|
| 774 | %/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 775 | global ZERO ONE |
---|
| 776 | global S4B sdMethod rkS rkGamma rkA rkB rkELO rkBhat |
---|
| 777 | global rkC rkD rkE rkTheta rkAlpha |
---|
| 778 | |
---|
| 779 | sdMethod = S4B; |
---|
| 780 | |
---|
| 781 | %/* Number of stages */ |
---|
| 782 | rkS = 5; |
---|
| 783 | |
---|
| 784 | %/* Method coefficients */ |
---|
| 785 | rkGamma = 0.25; |
---|
| 786 | |
---|
| 787 | rkA(1,1) = 0.25; |
---|
| 788 | rkA(1,2) = 0.5; |
---|
| 789 | rkA(2,2) = 0.25; |
---|
| 790 | rkA(1,3) = 0.34; |
---|
| 791 | rkA(2,3) = -0.40e-01; |
---|
| 792 | rkA(3,3) = 0.25; |
---|
| 793 | rkA(1,4) = 0.2727941176470588235294117647058824; |
---|
| 794 | rkA(2,4) = -0.5036764705882352941176470588235294e-01; |
---|
| 795 | rkA(3,4) = 0.2757352941176470588235294117647059e-01; |
---|
| 796 | rkA(4,4) = 0.25; |
---|
| 797 | rkA(1,5) = 1.041666666666666666666666666666667; |
---|
| 798 | rkA(2,5) = -1.020833333333333333333333333333333; |
---|
| 799 | rkA(3,5) = 7.812500000000000000000000000000000; |
---|
| 800 | rkA(4,5) = -7.083333333333333333333333333333333; |
---|
| 801 | rkA(5,5) = 0.25; |
---|
| 802 | |
---|
| 803 | rkB(1) = 1.041666666666666666666666666666667; |
---|
| 804 | rkB(2) = -1.020833333333333333333333333333333; |
---|
| 805 | rkB(3) = 7.812500000000000000000000000000000; |
---|
| 806 | rkB(4) = -7.083333333333333333333333333333333; |
---|
| 807 | rkB(5) = 0.250000000000000000000000000000000; |
---|
| 808 | |
---|
| 809 | rkBhat(1) = 1.069791666666666666666666666666667; |
---|
| 810 | rkBhat(2) = -0.894270833333333333333333333333333; |
---|
| 811 | rkBhat(3) = 7.695312500000000000000000000000000; |
---|
| 812 | rkBhat(4) = -7.083333333333333333333333333333333; |
---|
| 813 | rkBhat(5) = 0.212500000000000000000000000000000; |
---|
| 814 | |
---|
| 815 | rkC(1) = 0.25; |
---|
| 816 | rkC(2) = 0.75; |
---|
| 817 | rkC(3) = 0.55; |
---|
| 818 | rkC(4) = 0.5; |
---|
| 819 | rkC(5) = ONE; |
---|
| 820 | |
---|
| 821 | %/* Ynew = Yold + h*Sum_i {rkB_i*k_i} = Yold + Sum_i {rkD_i*Z_i} */ |
---|
| 822 | rkD(1) = ZERO; |
---|
| 823 | rkD(2) = ZERO; |
---|
| 824 | rkD(3) = ZERO; |
---|
| 825 | rkD(4) = ZERO; |
---|
| 826 | rkD(5) = ONE; |
---|
| 827 | |
---|
| 828 | %/* Err = h * Sum_i {(rkB_i-rkBhat_i)*k_i} = Sum_i {rkE_i*Z_i} */ |
---|
| 829 | rkE(1) = 0.5750; |
---|
| 830 | rkE(2) = 0.2125; |
---|
| 831 | rkE(3) = -4.6875; |
---|
| 832 | rkE(4) = 4.2500; |
---|
| 833 | rkE(5) = 0.1500; |
---|
| 834 | |
---|
| 835 | %/* Local order of Err estimate */ |
---|
| 836 | rkELO = 4; |
---|
| 837 | |
---|
| 838 | %/* h*Sum_j {rkA_ij*k_j} = Sum_j {rkTheta_ij*Z_j} */ |
---|
| 839 | rkTheta(1,2) = 2.0; |
---|
| 840 | rkTheta(1,3) = 1.680000000000000000000000000000000; |
---|
| 841 | rkTheta(2,3) = -0.1600000000000000000000000000000000; |
---|
| 842 | rkTheta(1,4) = 1.308823529411764705882352941176471; |
---|
| 843 | rkTheta(2,4) = -0.1838235294117647058823529411764706; |
---|
| 844 | rkTheta(3,4) = 0.1102941176470588235294117647058824; |
---|
| 845 | rkTheta(1,5) = -3.083333333333333333333333333333333; |
---|
| 846 | rkTheta(2,5) = -4.291666666666666666666666666666667; |
---|
| 847 | rkTheta(3,5) = 34.37500000000000000000000000000000; |
---|
| 848 | rkTheta(4,5) = -28.3333333333333333333333333333; |
---|
| 849 | |
---|
| 850 | %/* Starting value for Newton iterations: Z_i^0 = Sum_j {rkAlpha_ij*Z_j} */ |
---|
| 851 | rkAlpha(1,2) = 3.0; |
---|
| 852 | rkAlpha(1,3) = 0.8800000000000000000000000000000000; |
---|
| 853 | rkAlpha(2,3) = 0.4400000000000000000000000000000000; |
---|
| 854 | rkAlpha(1,4) = 0.1666666666666666666666666666666667; |
---|
| 855 | rkAlpha(2,4) = -0.8333333333333333333333333333333333e-01; |
---|
| 856 | rkAlpha(3,4) = 0.9469696969696969696969696969696970; |
---|
| 857 | rkAlpha(1,5) = -6.0; |
---|
| 858 | rkAlpha(2,5) = 9.0; |
---|
| 859 | rkAlpha(3,5) = -56.81818181818181818181818181818182; |
---|
| 860 | rkAlpha(4,5) = 54.0; |
---|
| 861 | |
---|
| 862 | %/* end Sdirk4b */ |
---|
| 863 | |
---|
| 864 | %/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 865 | function Sdirk2a() |
---|
| 866 | %/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 867 | global ZERO ONE |
---|
| 868 | global S2A sdMethod rkS rkGamma rkA rkB rkELO rkBhat |
---|
| 869 | global rkC rkD rkE rkTheta rkAlpha |
---|
| 870 | |
---|
| 871 | sdMethod = S2A; |
---|
| 872 | |
---|
| 873 | %/* ~~~> Number of stages */ |
---|
| 874 | rkS = 2; |
---|
| 875 | |
---|
| 876 | %/* ~~~> Method coefficients */ |
---|
| 877 | rkGamma = 0.2928932188134524755991556378951510; |
---|
| 878 | rkA(1,1) = 0.2928932188134524755991556378951510; |
---|
| 879 | rkA(1,2) = 0.7071067811865475244008443621048490; |
---|
| 880 | rkA(2,2) = 0.2928932188134524755991556378951510; |
---|
| 881 | rkB(1) = 0.7071067811865475244008443621048490; |
---|
| 882 | rkB(2) = 0.2928932188134524755991556378951510; |
---|
| 883 | rkBhat(1) = 0.6666666666666666666666666666666667; |
---|
| 884 | rkBhat(2) = 0.3333333333333333333333333333333333; |
---|
| 885 | rkC(1) = 0.292893218813452475599155637895151; |
---|
| 886 | rkC(2) = ONE; |
---|
| 887 | |
---|
| 888 | %/* ~~~> Ynew = Yold + h*Sum_i {rkB_i*k_i} = Yold + Sum_i {rkD_i*Z_i} */ |
---|
| 889 | rkD(1) = ZERO; |
---|
| 890 | rkD(2) = ONE; |
---|
| 891 | |
---|
| 892 | %/* ~~~> Err = h * Sum_i {(rkB_i-rkBhat_i)*k_i} = Sum_i {rkE_i*Z_i} */ |
---|
| 893 | rkE(1) = 0.4714045207910316829338962414032326; |
---|
| 894 | rkE(2) = -0.1380711874576983496005629080698993; |
---|
| 895 | |
---|
| 896 | %/* ~~~> Local order of Err estimate */ |
---|
| 897 | rkELO = 2; |
---|
| 898 | |
---|
| 899 | %/* ~~~> h*Sum_j {rkA_ij*k_j} = Sum_j {rkTheta_ij*Z_j} */ |
---|
| 900 | rkTheta(1,2) = 2.414213562373095048801688724209698; |
---|
| 901 | |
---|
| 902 | %/* ~~~> Starting value for Newton iterations */ |
---|
| 903 | rkAlpha(1,2) = 3.414213562373095048801688724209698; |
---|
| 904 | |
---|
| 905 | % /* end Sdirk2a */ |
---|
| 906 | |
---|
| 907 | %/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 908 | function Sdirk2b() |
---|
| 909 | %/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 910 | global ZERO ONE |
---|
| 911 | global S2B sdMethod rkS rkGamma rkA rkB rkELO rkBhat |
---|
| 912 | global rkC rkD rkE rkTheta rkAlpha |
---|
| 913 | |
---|
| 914 | sdMethod = S2B; |
---|
| 915 | |
---|
| 916 | %/* ~~~> Number of stages */ |
---|
| 917 | rkS = 2; |
---|
| 918 | |
---|
| 919 | %/* ~~~> Method coefficients */ |
---|
| 920 | rkGamma = 1.707106781186547524400844362104849; |
---|
| 921 | rkA(1,1) = 1.707106781186547524400844362104849; |
---|
| 922 | rkA(1,2) = -0.707106781186547524400844362104849; |
---|
| 923 | rkA(2,2) = 1.707106781186547524400844362104849; |
---|
| 924 | rkB(1) = -0.707106781186547524400844362104849; |
---|
| 925 | rkB(2) = 1.707106781186547524400844362104849; |
---|
| 926 | rkBhat(1) = 0.6666666666666666666666666666666667; |
---|
| 927 | rkBhat(2) = 0.3333333333333333333333333333333333; |
---|
| 928 | rkC(1) = 1.707106781186547524400844362104849; |
---|
| 929 | rkC(2) = ONE; |
---|
| 930 | |
---|
| 931 | %/* ~~~> Ynew = Yold + h*Sum_i {rkB_i*k_i} = Yold + Sum_i {rkD_i*Z_i} */ |
---|
| 932 | rkD(1) = ZERO; |
---|
| 933 | rkD(2) = ONE; |
---|
| 934 | |
---|
| 935 | %/* ~~~> Err = h * Sum_i {(rkB_i-rkBhat_i)*k_i} = Sum_i {rkE_i*Z_i} */ |
---|
| 936 | rkE(1) = -0.4714045207910316829338962414032326; |
---|
| 937 | rkE(2) = 0.8047378541243650162672295747365659; |
---|
| 938 | |
---|
| 939 | %/* ~~~> Local order of Err estimate */ |
---|
| 940 | rkELO = 2; |
---|
| 941 | |
---|
| 942 | %/* ~~~> h*Sum_j {rkA_ij*k_j} = Sum_j {rkTheta_ij*Z_j} */ |
---|
| 943 | rkTheta(1,2) = -0.414213562373095048801688724209698; |
---|
| 944 | |
---|
| 945 | %/* ~~~> Starting value for Newton iterations */ |
---|
| 946 | rkAlpha(1,2) = 0.5857864376269049511983112757903019; |
---|
| 947 | |
---|
| 948 | % /* end Sdirk2b */ |
---|
| 949 | |
---|
| 950 | %/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 951 | function Sdirk3a() |
---|
| 952 | %/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 953 | global ZERO ONE |
---|
| 954 | global S3A sdMethod rkS rkGamma rkA rkB rkELO rkBhat |
---|
| 955 | global rkC rkD rkE rkTheta rkAlpha |
---|
| 956 | |
---|
| 957 | sdMethod = S3A; |
---|
| 958 | |
---|
| 959 | %/* ~~~> Number of stages */ |
---|
| 960 | rkS = 3; |
---|
| 961 | |
---|
| 962 | %/* ~~~> Method coefficients */ |
---|
| 963 | rkGamma = 0.2113248654051871177454256097490213; |
---|
| 964 | rkA(1,1) = 0.2113248654051871177454256097490213; |
---|
| 965 | rkA(1,2) = 0.2113248654051871177454256097490213; |
---|
| 966 | rkA(2,2) = 0.2113248654051871177454256097490213; |
---|
| 967 | rkA(1,3) = 0.2113248654051871177454256097490213; |
---|
| 968 | rkA(2,3) = 0.5773502691896257645091487805019573; |
---|
| 969 | rkA(3,3) = 0.2113248654051871177454256097490213; |
---|
| 970 | rkB(1) = 0.2113248654051871177454256097490213; |
---|
| 971 | rkB(2) = 0.5773502691896257645091487805019573; |
---|
| 972 | rkB(3) = 0.2113248654051871177454256097490213; |
---|
| 973 | rkBhat(1)= 0.2113248654051871177454256097490213; |
---|
| 974 | rkBhat(2)= 0.6477918909913548037576239837516312; |
---|
| 975 | rkBhat(3)= 0.1408832436034580784969504064993475; |
---|
| 976 | rkC(1) = 0.2113248654051871177454256097490213; |
---|
| 977 | rkC(2) = 0.4226497308103742354908512194980427; |
---|
| 978 | rkC(3) = ONE; |
---|
| 979 | |
---|
| 980 | %/* ~~~> Ynew = Yold + h*Sum_i {rkB_i*k_i} = Yold + Sum_i {rkD_i*Z_i} */ |
---|
| 981 | rkD(1) = ZERO; |
---|
| 982 | rkD(2) = ZERO; |
---|
| 983 | rkD(3) = ONE; |
---|
| 984 | |
---|
| 985 | %/* ~~~> Err = h * Sum_i {(rkB_i-rkBhat_i)*k_i} = Sum_i {rkE_i*Z_i} */ |
---|
| 986 | rkE(1) = 0.9106836025229590978424821138352906; |
---|
| 987 | rkE(2) = -1.244016935856292431175815447168624; |
---|
| 988 | rkE(3) = 0.3333333333333333333333333333333333; |
---|
| 989 | |
---|
| 990 | %/* ~~~> Local order of Err estimate */ |
---|
| 991 | rkELO = 2; |
---|
| 992 | |
---|
| 993 | %/* ~~~> h*Sum_j {rkA_ij*k_j} = Sum_j {rkTheta_ij*Z_j} */ |
---|
| 994 | rkTheta(1,2) = ONE; |
---|
| 995 | rkTheta(1,3) = -1.732050807568877293527446341505872; |
---|
| 996 | rkTheta(2,3) = 2.732050807568877293527446341505872; |
---|
| 997 | |
---|
| 998 | %/* ~~~> Starting value for Newton iterations */ |
---|
| 999 | rkAlpha(1,2) = 2.0; |
---|
| 1000 | rkAlpha(1,3) = -12.92820323027550917410978536602349; |
---|
| 1001 | rkAlpha(2,3) = 8.83012701892219323381861585376468; |
---|
| 1002 | |
---|
| 1003 | % /* end Sdirk3a */ |
---|
| 1004 | |
---|
| 1005 | |
---|
| 1006 | |
---|