[2696] | 1 | function [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 | %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 108 | global Nfun Njac Nstp Nacc Nrej Ndec Nsol Nsng Ntexit Nhexit Nhnew |
---|
| 109 | % Parse ODE options |
---|
| 110 | AbsTol = odeget(Options, 'AbsTol'); |
---|
| 111 | if isempty(AbsTol) |
---|
| 112 | AbsTol = 1.0e-3; |
---|
| 113 | end |
---|
| 114 | Hstart = odeget(Options, 'InitialStep'); |
---|
| 115 | if isempty(Hstart) |
---|
| 116 | Hstart = 0; |
---|
| 117 | end |
---|
| 118 | Jacobian = odeget(Options, 'Jacobian'); |
---|
| 119 | if isempty(Jacobian) |
---|
| 120 | error('A Jacobian function is required.'); |
---|
| 121 | end |
---|
| 122 | Hmax = odeget(Options, 'MaxStep'); |
---|
| 123 | if isempty(Hmax) |
---|
| 124 | Hmax = 0; |
---|
| 125 | end |
---|
| 126 | RelTol = odeget(Options, 'RelTol'); |
---|
| 127 | if isempty(RelTol) |
---|
| 128 | RelTol = 1.0e-4; |
---|
| 129 | end |
---|
| 130 | % Initialize statistics |
---|
| 131 | Nfun=1; Njac=2; Nstp=3; Nacc=4; Nrej=5; Ndec=6; Nsol=7; Nsng=8; |
---|
| 132 | Ntexit=1; Nhexit=2; Nhnew=3; |
---|
| 133 | RSTATUS = zeros(20,1); |
---|
| 134 | ISTATUS = zeros(20,1); |
---|
| 135 | |
---|
| 136 | % Get problem size |
---|
| 137 | steps = length(Tspan); |
---|
| 138 | N = max(size(Y0)); |
---|
| 139 | |
---|
| 140 | % Initialize tolerances |
---|
| 141 | ATOL = ones(N,1)*AbsTol; |
---|
| 142 | RTOL = ones(N,1)*RelTol; |
---|
| 143 | |
---|
| 144 | % Initialize step |
---|
| 145 | RCNTRL(2) = max(0, Hmax); |
---|
| 146 | RCNTRL(3) = max(0, Hstart); |
---|
| 147 | |
---|
| 148 | % Integrate |
---|
| 149 | Y = zeros(N,steps); |
---|
| 150 | T = zeros(steps,1); |
---|
| 151 | Y(:,1) = Y0; |
---|
| 152 | T(1) = Tspan(1); |
---|
| 153 | t=cputime; |
---|
| 154 | for 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 | |
---|
| 163 | end |
---|
| 164 | cputime-t |
---|
| 165 | RSTAT=RSTATUS; |
---|
| 166 | ISTAT=ISTATUS; |
---|
| 167 | %~~~> Statistics on the work performed by the RK method |
---|
| 168 | function [T,Y,Ierr,ISTATUS,RSTATUS] = RK_Integrate(N,Y0,TIN, TOUT,OdeFunction1,... |
---|
| 169 | OdeJacobian1,ATOL,RTOL,RCNTRL,ICNTRL,ISTATUS,RSTATUS) |
---|
| 170 | global ZERO ONE Nfun Njac Nstp Nacc Nrej Ndec Nsol Nsng Ntexit Nhacc Nhnew |
---|
| 171 | global R2A R1A L3C GAU L3A RKmax |
---|
| 172 | ZERO =0;ONE=1;Nfun=1;Njac=2;Nstp=3;Nacc=4;Nrej=5;Ndec=6;Nsol=7;Nsng=8; |
---|
| 173 | Ntexit=1;Nhacc=2;Nhnew=3;RKmax=3;R2A=1;R1A=2;L3C=3;GAU=4;L3A=5; |
---|
| 174 | Ierr =0; |
---|
| 175 | RSTATUS = zeros(20,1); |
---|
| 176 | ISTATUS = zeros(20,1); |
---|
| 177 | T1 = TIN; |
---|
| 178 | T2 = TOUT; |
---|
| 179 | Y(:,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 | |
---|
| 186 | if(Ierr < 0) |
---|
| 187 | disp('Runge-kutta Unsuccessful exit at T=',num2str(TIN)); |
---|
| 188 | end |
---|
| 189 | return |
---|
| 190 | |
---|
| 191 | function [T,Y,Ierr,RSTATUS,ISTATUS]=RungeKutta(N,T1,T2,Y,RelTol,AbsTol,RCNTRL,... |
---|
| 192 | ICNTRL,ISTATUS,RSTATUS,Ierr,OdeFunction1,OdeJacobian1) |
---|
| 193 | global ZERO Roundoff SdirkError |
---|
| 194 | T=T1; |
---|
| 195 | |
---|
| 196 | if (ICNTRL(1)==0) |
---|
| 197 | ITOL=1; |
---|
| 198 | else |
---|
| 199 | ITOL=0; |
---|
| 200 | end |
---|
| 201 | |
---|
| 202 | if(ICNTRL(9)==0) |
---|
| 203 | SdirkError = 0; |
---|
| 204 | else |
---|
| 205 | SdirkError = 1; |
---|
| 206 | end |
---|
| 207 | |
---|
| 208 | switch(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; |
---|
| 223 | end |
---|
| 224 | |
---|
| 225 | if (ICNTRL(3)==0) |
---|
| 226 | Max_no_steps=200000; |
---|
| 227 | else |
---|
| 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 |
---|
| 232 | end |
---|
| 233 | |
---|
| 234 | if (ICNTRL(4)==0) |
---|
| 235 | NewtonMaxit=8; |
---|
| 236 | else |
---|
| 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 |
---|
| 242 | end |
---|
| 243 | |
---|
| 244 | if (ICNTRL(5)==0) |
---|
| 245 | StartNewton=1; |
---|
| 246 | else |
---|
| 247 | StartNewton =0; |
---|
| 248 | end |
---|
| 249 | |
---|
| 250 | if(ICNTRL(10)==0) |
---|
| 251 | Gustafsson = 1; |
---|
| 252 | else |
---|
| 253 | Gustafsson = 0; |
---|
| 254 | end |
---|
| 255 | Roundoff = eps; |
---|
| 256 | if(RCNTRL(1)==ZERO) |
---|
| 257 | Hmin = ZERO; |
---|
| 258 | else |
---|
| 259 | Hmin = min(abs(RCNTRL(1)),abs(T2 -T1)); |
---|
| 260 | end |
---|
| 261 | if(RCNTRL(2)==ZERO) |
---|
| 262 | Hmax=abs(T2-T1); |
---|
| 263 | else |
---|
| 264 | Hmax= min(abs(RCNTRL(2)),abs(T2-T1)); |
---|
| 265 | end |
---|
| 266 | if RCNTRL(3)==0 |
---|
| 267 | Hstart=0; |
---|
| 268 | else |
---|
| 269 | Hstart = min(abs(RCNTRL(3)),abs(T2-T1)); |
---|
| 270 | end |
---|
| 271 | |
---|
| 272 | %FacMin-- Lower bound on step decrease factor |
---|
| 273 | if (RCNTRL(4) == 0) |
---|
| 274 | FacMin = 0.2; |
---|
| 275 | else |
---|
| 276 | FacMin = RCNTRL(4); |
---|
| 277 | end |
---|
| 278 | %FacMax--Upper bound on step increase factor |
---|
| 279 | if RCNTRL(5)==0 |
---|
| 280 | FacMax=8.0; |
---|
| 281 | else |
---|
| 282 | FacMax=RCNTRL(5); |
---|
| 283 | end |
---|
| 284 | %FacRej--step decrease factor after 2 consecutive rejections |
---|
| 285 | if RCNTRL(6)==0 |
---|
| 286 | FacRej=0.1; |
---|
| 287 | else |
---|
| 288 | FacRej=RCNTRL(6); |
---|
| 289 | end |
---|
| 290 | |
---|
| 291 | %Facsafe:by which the new step is slightly smaller than the |
---|
| 292 | %predicted value |
---|
| 293 | if RCNTRL(7)==0 |
---|
| 294 | FacSafe=0.9; |
---|
| 295 | else |
---|
| 296 | FacSafe=RCNTRL(7); |
---|
| 297 | end |
---|
| 298 | |
---|
| 299 | if ( (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); |
---|
| 302 | end |
---|
| 303 | |
---|
| 304 | %ThetaMin: Decides whether the jacobian should be recomputed |
---|
| 305 | if RCNTRL(8)==0 |
---|
| 306 | ThetaMin=1.0e-03; |
---|
| 307 | else |
---|
| 308 | ThetaMin=RCNTRL(8); |
---|
| 309 | end |
---|
| 310 | if (ThetaMin <= 0.0 || ThetaMin >= 1.0) |
---|
| 311 | disp(['RCNTRL[8]=', num2str(RCNTRL(8))]); |
---|
| 312 | RK_ErrorMsg(-5, Tin, ZERO, Ierr); |
---|
| 313 | end |
---|
| 314 | |
---|
| 315 | if RCNTRL(9)==0 |
---|
| 316 | NewtonTol = 3.0e-02; |
---|
| 317 | else |
---|
| 318 | NewtonTol = RCNTRL(9); |
---|
| 319 | if NewtonTol<=Roundoff |
---|
| 320 | disp(['RCNTRL(9)=',num2str(RCNTRL(9))]); |
---|
| 321 | RK_ErrorMsg(-6, Tin, ZERO,Ierr); |
---|
| 322 | end |
---|
| 323 | end |
---|
| 324 | |
---|
| 325 | %Qmin and Qmax: If Qmin < Hnew/Hold < Qmax then step size is constant |
---|
| 326 | if RCNTRL(10)==0 |
---|
| 327 | Qmin=1; |
---|
| 328 | else |
---|
| 329 | Qmin = RCNTRL(10); |
---|
| 330 | end |
---|
| 331 | if RCNTRL(11)==0 |
---|
| 332 | Qmax=1.2; |
---|
| 333 | else |
---|
| 334 | Qmax=RCNTRL(11); |
---|
| 335 | end |
---|
| 336 | if (Qmin > 1 || Qmax <1) |
---|
| 337 | disp(['Qmin=', num2str(RCNTRL(10)) ' Qmax=',num2str(RCNTRl(11))]); |
---|
| 338 | RK_ErrorMsg(-7, T, ZERO, Ierr); |
---|
| 339 | end |
---|
| 340 | |
---|
| 341 | %check if tolerances are reasonable |
---|
| 342 | if 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 |
---|
| 347 | else |
---|
| 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 |
---|
| 355 | end |
---|
| 356 | |
---|
| 357 | if(Ierr < 0) |
---|
| 358 | return; |
---|
| 359 | end |
---|
| 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); |
---|
| 365 | return; |
---|
| 366 | |
---|
| 367 | |
---|
| 368 | |
---|
| 369 | function [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) |
---|
| 374 | global SdirkError ZERO Nfun Njac Nsng Nstp ONE rkMethod L3A rkBgam |
---|
| 375 | global rkTheta rkGamma rkD rkELO Nacc Ntexit Nhacc Nhnew Hacc |
---|
| 376 | T = T1; |
---|
| 377 | CONT=zeros(N,3); |
---|
| 378 | Tdirection = sign(T2-T); |
---|
| 379 | H = min (max(abs(Hmin),abs(Hstart)),Hmax); |
---|
| 380 | if(abs(H)<=10*Roundoff) |
---|
| 381 | H=1.0e-6; |
---|
| 382 | end |
---|
| 383 | H =sign(Tdirection)*abs(H); |
---|
| 384 | Hold = H; |
---|
| 385 | Reject = 0; |
---|
| 386 | FirstStep = 1; |
---|
| 387 | SkipJac = 0; |
---|
| 388 | SkipLU = 0; |
---|
| 389 | if((T+H*1.0001-T2)*Tdirection>=ZERO) |
---|
| 390 | H = T2 -T; |
---|
| 391 | end |
---|
| 392 | Nconsecutive = 0; |
---|
| 393 | SCAL = RK_ErrorScale(N,ITOL,AbsTol,RelTol,Y); |
---|
| 394 | while((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 | |
---|
| 661 | return; |
---|
| 662 | |
---|
| 663 | function [Err,ISTATUS]=RK_ErrorEstimate(N,H,T,Y,F0,E1L,E1U,Z1,Z2,Z3,SCAL,... |
---|
| 664 | FirstStep,Reject,ISTATUS) |
---|
| 665 | global rkE rkMethod R1A GAU L3A |
---|
| 666 | HrkE1=rkE(2)/H; |
---|
| 667 | HrkE2=rkE(3)/H; |
---|
| 668 | HrkE3=rkE(4)/H; |
---|
| 669 | F2=HrkE1*Z1+HrkE2*Z2+HrkE3*Z3; |
---|
| 670 | TMP = rkE(1)*F0+F2; |
---|
| 671 | [TMP,ISTATUS]=RK_KppSolve(E1L,E1U,TMP,ISTATUS); |
---|
| 672 | if((rkMethod == R1A) || (rkMethod == GAU) || (rkMethod == L3A)) |
---|
| 673 | [TMP,ISTATUS]=RK_KppSolve(E1L,E1U,TMP,ISTATUS); |
---|
| 674 | end |
---|
| 675 | |
---|
| 676 | if (rkMethod == GAU) |
---|
| 677 | [TMP,ISTATUS]=RK_KppSolve(E1L,E1U,TMP,ISTATUS); |
---|
| 678 | end |
---|
| 679 | Err=RK_ErrorNorm(N,SCAL,TMP); |
---|
| 680 | if(Err < 1) |
---|
| 681 | return; |
---|
| 682 | end |
---|
| 683 | if((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); |
---|
| 690 | end |
---|
| 691 | return |
---|
| 692 | |
---|
| 693 | function 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)]); |
---|
| 729 | return; |
---|
| 730 | |
---|
| 731 | function 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; |
---|
| 737 | return; |
---|
| 738 | |
---|
| 739 | function ErrorNorm = RK_ErrorNorm(N,SCAL,DY) |
---|
| 740 | |
---|
| 741 | temp=DY.*DY.*SCAL.*SCAL; |
---|
| 742 | ErrorNorm = sum(temp); |
---|
| 743 | ErrorNorm=max(sqrt(ErrorNorm/N),1.0e-10); |
---|
| 744 | |
---|
| 745 | function [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; |
---|
| 780 | return; |
---|
| 781 | |
---|
| 782 | |
---|
| 783 | |
---|
| 784 | function [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; |
---|
| 813 | return; |
---|
| 814 | |
---|
| 815 | function [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 | |
---|
| 848 | return |
---|
| 849 | |
---|
| 850 | function [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; |
---|
| 879 | return |
---|
| 880 | |
---|
| 881 | function [R4,ISTATUS]=RK_KppSolve(E1L,E1U,R4,ISTATUS) |
---|
| 882 | global Nsol |
---|
| 883 | temp = E1L\R4; |
---|
| 884 | R4=E1U\temp; |
---|
| 885 | ISTATUS(Nsol)=ISTATUS(Nsol)+1; |
---|
| 886 | return |
---|
| 887 | |
---|
| 888 | |
---|
| 889 | |
---|
| 890 | |
---|
| 891 | function 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; |
---|
| 1011 | return; |
---|
| 1012 | |
---|
| 1013 | function 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; |
---|
| 1133 | end; |
---|
| 1134 | return; |
---|
| 1135 | |
---|
| 1136 | % /* Lobatto3C_Coefficients */ |
---|
| 1137 | function 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; |
---|
| 1255 | return; |
---|
| 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 | %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1264 | function 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; |
---|
| 1378 | return; |
---|
| 1379 | % /* Radau1A_Coefficients */ |
---|
| 1380 | |
---|
| 1381 | %/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 1382 | % The coefficients of the 4-stage Lobatto-3A method |
---|
| 1383 | % (given to ~30 accurate digits) |
---|
| 1384 | % ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 1385 | function 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; |
---|
| 1517 | return; |
---|
| 1518 | %/* Lobatto3A_Coefficients */ |
---|
| 1519 | |
---|