1function varargout = ode2345(svname,ode,tspan,y0,options,varargin)
2
3
4
5
6
7
8
9GraphicalOutput = nargout == 0;
10StructOutput = nargout == 1;
11ArrayOutput = ~GraphicalOutput && ~StructOutput;
12
13
14eml_invariant(ArrayOutput, ...
15 'Coder:toolbox:RequiresNOutputs',svname,2);
16
17
18
19
20eml_invariant(eml_option('VariableSizing'), ...
21 'Coder:toolbox:RequiresVariableSizing',svname);
22
23
24
25eml_invariant(nargin >= 5,'MATLAB:minrhs');
26coder.internal.prefer_const(svname,ode,tspan,options);
27eml_invariant(eml_is_const(svname) && ...
28 (strcmp(svname,'ode23') || strcmp(svname,'ode45')), ...
29 'Coder:toolbox:UnsupportedSolver');
30ODE23 = coder.const(strcmp(svname,'ode23'));
31eml_invariant((ischar(ode) && eml_is_const(ode)) || ...
32 isa(ode,'function_handle'), ...
33 'Coder:toolbox:UnsupportedODEFunction');
34if ischar(ode)
35 ODEFunction = str2func(ode);
36else
37 ODEFunction = ode;
38end
39eml_invariant(isreal(tspan), ...
40 'Coder:toolbox:TspanMustBeReal');
41coder.internal.errorIf(isempty(tspan) || isempty(y0), ...
42 'MATLAB:odearguments:TspanOrY0NotSupplied',svname);
43coder.internal.errorIf(length(tspan) < 2, ...
44 'MATLAB:odearguments:SizeTspan',svname);
45t0 = tspan(1);
46tfinal = tspan(end);
47coder.internal.errorIf(t0 == tfinal, ...
48 'MATLAB:odearguments:TspanEndpointsNotDistinct');
49eml_invariant(ismonotonic(tspan), ...
50 'MATLAB:odearguments:TspanNotMonotonic');
51eml_invariant(isvector(y0), ...
52 'Coder:toolbox:y0MustBeAVector');
53
54
55
56Moption = odeget(options,'Mass','','fast');
57if ~isempty(Moption)
58 Msingular = odeget(options,'MassSingular','no','fast');
59 coder.internal.errorIf(strcmp(Msingular,'yes'), ...
60 make_error_id(svname,'MassSingularYes'));
61 if strcmp(Msingular,'maybe')
62 coder.internal.warning(make_error_id(svname,'MassSingularAssumedNo'));
63 end
64end
65
66
67
68ZERO = coder.internal.indexInt(0);
69ONE = coder.internal.indexInt(1);
70nsteps = ZERO;
71nfailed = ZERO;
72nfevals = ZERO;
73neq = coder.internal.indexInt(numel(y0));
74
75
76
77
78idxNonNegative = getIdxNonNegative(options,y0);
79AnyNonNegative = ~isempty(idxNonNegative);
80
81
82eml_invariant(isfloat(y0) && isa(y0,class(t0)), ...
83 'Coder:toolbox:InconsistentTypes');
84f0 = callODEFunctionNSM(ODEFunction,t0,y0(:),options,varargin{:});
85nfevals = nfevals + 1;
86datatype = coder.internal.scalarEg(tspan,y0,f0);
87realtype = real(datatype);
88realnull = zeros(0,0,'like',realtype);
89datanull = zeros(0,0,'like',datatype);
90outcls = class(datatype);
91
92
93eml_invariant(isfloat(datatype) && ...
94 isa(tspan,outcls) && isa(y0,outcls) && isa(f0,outcls), ...
95 'Coder:toolbox:InconsistentTypes');
96[m,n] = size(f0);
97
98eml_invariant(eml_is_const(n) && n == 1, ...
99 'Coder:toolbox:FOMustReturnCol');
100
101eml_invariant(m == neq, ...
102 'Coder:toolbox:SizeIC')
103
104
105
106rtol = cast(odeget(options,'RelTol',1e-3,'fast'),'like',realtype);
107eml_invariant(isscalar(rtol) && rtol > 0, ...
108 'MATLAB:odearguments:RelTolNotPosScalar');
109minrtol = 100*eps(outcls);
110if rtol < minrtol
111 rtol = minrtol;
112 coder.internal.warning('MATLAB:odearguments:RelTolIncrease', ...
113 coder.internal.num2str(rtol));
114end
115atol = cast(odeget(options,'AbsTol',1e-6,'fast'),'like',realtype);
116eml_invariant(ispositive(atol), ...
117 'MATLAB:odearguments:AbsTolNotPos');
118normcontrol = strcmp(odeget(options,'NormControl','off','fast'),'on');
119if normcontrol
120 eml_invariant(isscalar(atol), ...
121 'MATLAB:odearguments:NonScalarAbsTol');
122 normy = norm(y0(:));
123else
124
125 eml_invariant(isscalar(atol) || numel(atol) == neq, ...
126 'Coder:toolbox:SizeAbsTol',neq);
127 normy = zeros('like',realtype);
128end
129threshold = atol(:)/rtol;
130
131
132
133if ODE23
134 nstages = coder.internal.indexInt(4);
135 pow = cast(1/3,'like',realtype);
136 A = cast([1/2, 3/4, 1],'like',realtype);
137
138
139
140
141
142
143
144 B = cast( ...
145 [1/2; ...
146 0; 3/4; ...
147 2/9; 1/3; 4/9], ...
148 outcls);
149 E = cast([-5/72; 1/12; 1/9; -1/8],'like',realtype);
150 ntrpfun = @ntrp23;
151 MinStepFactor = cast(0.5,'like',realtype);
152 DefaultRefine = 1;
153else
154 nstages = coder.internal.indexInt(7);
155 pow = cast(1/5,'like',realtype);
156 A = cast([1/5,3/10,4/5,8/9,1,1],'like',realtype);
157
158
159
160
161
162
163
164
165
166
167 B = cast( ...
168 [1/5; ...
169 3/40; 9/40; ...
170 44/45; -56/15; 32/9; ...
171 19372/6561; -25360/2187; 64448/6561; -212/729; ...
172 9017/3168; -355/33; 46732/5247; 49/176; -5103/18656; ...
173 35/384; 0; 500/1113; 125/192; -2187/6784; 11/84], ...
174 'like',realtype);
175 E = cast( ...
176 [71/57600; 0; -71/16695; 71/1920; -17253/339200; 22/525; -1/40],...
177 'like',realtype);
178 ntrpfun = @ntrp45;
179 MinStepFactor = cast(0.1,'like',realtype);
180 DefaultRefine = 4;
181end
182
183
184
185eventOption = odeget(options,'Events','','fast');
186HasEvents = ~isempty(eventOption);
187if HasEvents
188
189 eml_invariant(~strcmpi(eml_option('UseMalloc'),'Off'), ...
190 'Coder:toolbox:EventsRequireDynamicMemoryAllocation',svname);
191 if ischar(eventOption)
192 EventFunction = str2func(eventOption);
193 else
194 EventFunction = eventOption;
195 end
196 valt = EventFunction(t0,y0(:),varargin{:});
197end
198
199
200
201
202teout = zeros(1,0,'like',realtype);
203yeout = zeros(neq,0,'like',datatype);
204ieout = zeros(1,0,coder.internal.indexIntClass);
205ntspan = numel(tspan);
206refine = coder.internal.indexInt( ...
207 max(1,odeget(options,'Refine',DefaultRefine,'fast')));
208if eml_is_const(ntspan) && ntspan == 2
209 coder.varsize('tout',[],[0,1]);
210 coder.varsize('yout');
211end
212if eml_is_const(size(tspan)) && ntspan == 2
213
214 chunk = coder.internal.indexInt( ...
215 min(max(100,50*refine),refine+floor((2^13)/neq)));
216 tout = zeros(1,chunk,'like',realtype);
217 yout = zeros(neq,chunk,'like',datatype);
218 if refine <= 1
219 outputAt = SOLVERSTEPS;
220 else
221 outputAt = REFINEDSTEPS;
222 S = cast(1:refine-1,'like',realtype)/cast(refine,'like',realtype);
223 end
224else
225 outputAt = REQUESTEDPOINTS;
226 if ntspan == 2
227
228
229 coder.internal.warning('Coder:toolbox:VarsizeTspanHasLengthTwo');
230 end
231 tout = zeros(1,ntspan,'like',realtype);
232 yout = zeros(neq,ntspan,'like',datatype);
233end
234nout = ONE;
235tout(nout) = t0;
236yout(:,nout) = y0(:);
237
238
239
240outputFunctionOption = odeget(options,'OutputFcn','','fast');
241HasOutputFunction = ~isempty(outputFunctionOption);
242if HasOutputFunction
243 if ischar(outputFunctionOption)
244 OutputFunction = str2func(outputFunctionOption);
245 else
246 OutputFunction = outputFunctionOption;
247 end
248 outputidx = odeget(options,'OutputSel',ONE:neq,'fast');
249 eml_invariant(isValidIdxVector(outputidx,neq), ...
250 'Coder:toolbox:InvalidOutputSel');
251 noutputidx = coder.internal.indexInt(numel(outputidx));
252 OutputFunction([t0,tfinal],y0(outputidx(:)),'init',varargin{:});
253end
254
255
256
257PrintStats = strcmp(odeget(options,'Stats','off','fast'),'on');
258
259eml_invariant(~PrintStats, ...
260 'Coder:toolbox:NoPrintStats');
261
262
263
264
265hmax = min( ...
266 abs(tfinal - t0), ...
267 abs(odeget(options,'MaxStep',0.1*(tfinal - t0),'fast')));
268eml_invariant(hmax > 0, ...
269 'MATLAB:odearguments:MaxStepLEzero');
270htry = abs(odeget(options,'InitialStep',[],'fast'));
271eml_invariant(isempty(htry) || htry > 0, ...
272 'MATLAB:odearguments:InitialStepLEzero');
273hmin = 16*eps(t0);
274if isempty(htry)
275
276 htspan = abs(tspan(2) - tspan(1));
277 absh = min(hmax,htspan);
278 if normcontrol
279 rh = (norm(f0)/max(normy,threshold))/(0.8*rtol^pow);
280 else
281 rh = norm(f0./max(abs(y0(:)),threshold),'inf')/(0.8*rtol^pow);
282 end
283 if absh*rh > 1
284 absh = 1/rh;
285 end
286 absh = max(absh,hmin);
287else
288 absh = min(hmax,max(hmin,htry));
289end
290
291
292
293t = t0;
294y = cast(y0(:),'like',datatype);
295f = zeros(neq,nstages,'like',y);
296f(:,1) = f0;
297tdir = sign(tfinal - t0);
298next = coder.internal.indexInt(2);
299TerminalEvent = false;
300
301
302
303MinStepExit = false;
304Done = false;
305while ~Done
306
307
308 hmin = 16*eps(t);
309 absh = min(hmax,max(hmin,absh));
310 h = tdir*absh;
311
312 if 1.1*absh >= abs(tfinal - t)
313 h = tfinal - t;
314 absh = abs(h);
315 Done = true;
316 end
317
318 NoFailedAttempts = true;
319 NNresetLastStage = false;
320 while true
321 Bcolidx = ONE;
322 for j = 1:nstages-2
323
324 Bcolidx = Bcolidx + j - 1;
325 ystage = coder.internal.blas.xgemv( ...
326 'N',neq,j,h,f,1,neq,B,Bcolidx,1,1,y,1,1);
327 f(:,j+1) = callODEFunctionNSM( ...
328 ODEFunction,t+h*A(j),ystage,options,varargin{:});
329 end
330 tnew = t + h*A(nstages-1);
331 if Done
332 tnew = tfinal;
333 end
334
335 Bcolidx = Bcolidx + nstages - 2;
336 ynew = coder.internal.blas.xgemv( ...
337 'N',neq,nstages-1,h,f,1,neq,B,Bcolidx,1,1,y,1,1);
338 f(:,nstages) = callODEFunctionNSM( ...
339 ODEFunction,tnew,ynew,options,varargin{:});
340 nfevals = nfevals + (nstages - 1);
341 h = tnew - t;
342
343 if normcontrol
344 normynew = norm(ynew);
345 errwt = max3(normy,normynew,threshold);
346 err = absh*(norm(f*E)/errwt);
347 else
348 err = absh*norm((f*E)./maxAbsThresh(y,ynew,threshold),'inf');
349 errwt = threshold;
350 end
351 NNrejectStep = false;
352 if AnyNonNegative && err <= rtol
353 errNN = normNegative(ynew,errwt,idxNonNegative,normcontrol);
354 if errNN > rtol
355 err = errNN;
356 NNrejectStep = true;
357 end
358 end
359
360
361
362
363 if err > rtol
364 nfailed = nfailed + 1;
365 if absh <= hmin
366 coder.internal.warning( ...
367 make_error_id(svname,'IntegrationTolNotMet'), ...
368 coder.internal.num2str(t), ...
369 coder.internal.num2str(hmin));
370 MinStepExit = true;
371 break
372 else
373 if NoFailedAttempts
374 NoFailedAttempts = false;
375 if NNrejectStep
376 absh = max(hmin,0.5*absh);
377 else
378 absh = max(hmin, ...
379 absh*max(MinStepFactor,0.8*(rtol/err)^pow));
380 end
381 else
382 absh = max(hmin,0.5*absh);
383 end
384 h = tdir*absh;
385 Done = false;
386 end
387 else
388 NNresetLastStage = false;
389 if AnyNonNegative
390 for k = 1:numel(idxNonNegative)
391 if ynew(idxNonNegative(k)) < 0
392 ynew(idxNonNegative(k)) = 0;
393 NNresetLastStage = true;
394 end
395 end
396 if normcontrol
397 normynew = norm(ynew);
398 end
399 end
400 break
401 end
402 end
403 if MinStepExit
404 break
405 end
406 nsteps = nsteps + 1;
407 if HasEvents
408 [ne,te,ye,ie,valt,TerminalEvent(1)] = ...
409 odezero(ntrpfun,EventFunction,valt,t,y,tnew,ynew,t0,h,f, ...
410 idxNonNegative,varargin{:});
411 if ne > 0
412 if StructOutput || nargout > 2
413 teout = [teout,te];
414 yeout = [yeout,ye];
415 ieout = [ieout,double(ie)];
416 end
417 if TerminalEvent
418
419
420 taux = t + (te(end) - t)*A;
421 [~,f(:,2:nstages)] = ntrpfun(taux,t,y, ...
422 realnull,datanull,h,f,idxNonNegative);
423 tnew = te(end);
424 ynew = ye(:,end);
425 h = tnew - t;
426 Done = true;
427 end
428 end
429 end
430 switch outputAt
431 case SOLVERSTEPS
432 nout = nout + 1;
433 if nout > length(tout)
434 tout = [tout,zeros(1,chunk,'like',t)];
435 yout = [yout,zeros(neq,chunk,'like',y)];
436 end
437 tout(nout) = tnew;
438 for j = 1:neq
439 yout(j,nout) = ynew(j);
440 end
441 noutnew = ONE;
442 outidx = nout;
443 case REFINEDSTEPS
444 noutnew = refine;
445 outidx = nout + 1;
446 tref = t + (tnew - t)*S;
447 toutnew = [tref,tnew];
448 youtnew = [ntrpfun(tref,t,y, ...
449 realnull,datanull,h,f,idxNonNegative),ynew];
450 nout = nout + noutnew;
451 if nout > length(tout)
452
453 tout = [tout,zeros(1,chunk,'like',t)];
454 yout = [yout,zeros(neq,chunk,'like',y)];
455 end
456
457
458 for k = 1:noutnew
459 tout(k - 1 + outidx) = toutnew(k);
460 for j = 1:neq
461 yout(j,k - 1 + outidx) = youtnew(j,k);
462 end
463 end
464 case REQUESTEDPOINTS
465 nnxt = NextNext(next,ntspan,tspan,tdir,tnew);
466 noutnew = nnxt - next;
467 outidx = next;
468 tnewused = false;
469 if noutnew > 0
470 i2 = nnxt - 1;
471 for k = outidx:i2-1
472 tout(k) = tspan(k);
473 yout(:,k) = ntrpfun(tspan(k),t,y, ...
474 realnull,datanull,h,f,idxNonNegative);
475 end
476 tout(i2) = tspan(i2);
477 if tspan(i2) == tnew
478 yout(:,i2) = ynew;
479 tnewused = true;
480 else
481 yout(:,i2) = ntrpfun(tspan(i2),t,y, ...
482 realnull,datanull,h,f,idxNonNegative);
483 end
484 nout = nout + noutnew;
485 next = nnxt;
486 end
487 if HasEvents && TerminalEvent && ~tnewused && next <= ntspan
488
489 tout(next) = tnew;
490 yout(:,next) = ynew;
491 nout = nout + 1;
492 end
493 end
494 if HasOutputFunction
495 toutnew = zeros(1,noutnew,'like',realtype);
496 youtnew = zeros(noutputidx,noutnew,'like',datatype);
497 for k = 1:noutnew
498 toutnew(k) = tout(k - 1 + outidx);
499 for j = 1:noutputidx
500 youtnew(j,k) = yout(outputidx(j),k - 1 + outidx);
501 end
502 end
503 TerminalEvent(1) = OutputFunction( ...
504 toutnew,youtnew,'',varargin{:});
505 if TerminalEvent
506 Done = true;
507 end
508 end
509 if Done
510 break
511 end
512
513 if NoFailedAttempts
514
515 temp = 1.25*(err/rtol)^pow;
516 if temp > 0.2
517 absh = absh/temp;
518 else
519 absh = 5*absh;
520 end
521 end
522
523 t = tnew;
524 y = ynew;
525 if normcontrol
526 normy = normynew;
527 end
528 if NNresetLastStage
529
530
531 f(:,nstages) = callODEFunctionNSM( ...
532 ODEFunction,tnew,ynew,options,varargin{:});
533 nfevals = nfevals + 1;
534 end
535 f(:,1) = f(:,nstages);
536end
537eml_invariant(nargout <= 3 || (HasEvents && nargout <= 6), ...
538 'MATLAB:maxlhs');
539varargout{1} = tout(1:nout).';
540varargout{2} = yout(:,1:nout).';
541if HasEvents
542 varargout{3} = teout.';
543 varargout{4} = yeout.';
544 varargout{5} = double(ieout.');
545 varargout{6} = double([nsteps;nfailed;nfevals;0;0;0]);
546else
547 varargout{3} = double([nsteps;nfailed;nfevals;0;0;0]);
548
549 varargout{4} = [];
550 varargout{5} = [];
551 varargout{6} = [];
552end
553if HasOutputFunction
554 OutputFunction(realnull,datanull,'done',varargin{:});
555end