Mercurial > octave

annotate scripts/ode/private/runge_kutta_interpolate.m @ 20903:3d3da166dac5

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
2015 Code Sprint: finish import of ode23 into core * scripts/ode/private/runge_kutta_23.m: apply vectorization. * scripts/ode/private/runge_kutta_45_dorpri.m: remove unused parts of the tableau. * scripts/ode/private/runge_kutta_interpolate.m: reimplement cubic hermite interpolation.
author Carlo de Falco <carlo.defalco@polimi.it>
date Tue, 15 Dec 2015 18:50:58 +0100
parents ddc18b909ec7
children ebe061d6feea
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
20565
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
1 ## Copyright (C) 2015 Carlo de Falco
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
2 ## Copyright (C) 2015 Jacopo Corno <jacopo.corno@gmail.com>
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
3 ##
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
4 ## This file is part of Octave.
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
5 ##
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
6 ## Octave is free software; you can redistribute it and/or modify it
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
7 ## under the terms of the GNU General Public License as published by
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
8 ## the Free Software Foundation; either version 3 of the License, or (at
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
9 ## your option) any later version.
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
10 ##
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
11 ## Octave is distributed in the hope that it will be useful, but
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
12 ## WITHOUT ANY WARRANTY; without even the implied warranty of
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
14 ## General Public License for more details.
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
15 ##
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
16 ## You should have received a copy of the GNU General Public License
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
17 ## along with Octave; see the file COPYING. If not, see
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
18 ## <http://www.gnu.org/licenses/>.
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
19
20847
ddc18b909ec7 codesprint: ode: first and second order interpolator for dense_output
jcorno <jacopo.corno@gmail.com>
parents: 20634
diff changeset
20 function u_interp = runge_kutta_interpolate (order, z, u, t, k_vals, dt, func, args)
20565
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
21
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20567
diff changeset
22 switch (order)
20565
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
23
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20567
diff changeset
24 ## FIXME: Can interpolations for orders 1-4 simply be deleted? 2015-10-14.
20847
ddc18b909ec7 codesprint: ode: first and second order interpolator for dense_output
jcorno <jacopo.corno@gmail.com>
parents: 20634
diff changeset
25
20565
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
26 case 1
20847
ddc18b909ec7 codesprint: ode: first and second order interpolator for dense_output
jcorno <jacopo.corno@gmail.com>
parents: 20634
diff changeset
27 u_interp = interp1 (z, u', t, 'linear');
ddc18b909ec7 codesprint: ode: first and second order interpolator for dense_output
jcorno <jacopo.corno@gmail.com>
parents: 20634
diff changeset
28
20565
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
29 case 2
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
30 if (! isempty (k_vals))
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
31 der = k_vals(:,1);
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
32 else
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
33 der = feval (func, z(1) , u(:,1), args);
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
34 endif
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
35 u_interp = quadratic_interpolation (z, u, der, t);
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
36 case 3
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
37 u_interp = ...
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
38 hermite_cubic_interpolation (z, u, k_vals, t);
20903
3d3da166dac5 2015 Code Sprint: finish import of ode23 into core
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20847
diff changeset
39 #{
20565
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
40 case 4
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
41 ## if ode45 is used without local extrapolation this function
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
42 ## doesn't require a new function evaluation.
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
43 u_interp = dorpri_interpolation ([z(i-1) z(i)],
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
44 [u(:,i-1) u(:,i)],
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
45 k_vals, tspan(counter));
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
46 #}
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
47 case 5
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
48 ## ode45 with Dormand-Prince scheme:
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20567
diff changeset
49 u_interp = hermite_quartic_interpolation (z, u, k_vals, t);
20565
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
50
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
51 ## it is also possible to do a new function evaluation and use
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
52 ## the quintic hermite interpolator
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
53 ## f_half = feval (func, t+1/2*dt, u_half,
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
54 ## options.funarguments{:});
20565
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
55 ## u_interp =
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
56 ## hermite_quintic_interpolation ([z(i-1) z(i)],
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
57 ## [u(:,i-1) u_half u(:,i)],
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
58 ## [k_vals(:,1) f_half ...
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
59 ## k_vals(:,end)],
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
60 ## tspan(counter));
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20567
diff changeset
61
20565
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
62 otherwise
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20567
diff changeset
63 warning (["High order interpolation not yet implemented: ", ...
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20567
diff changeset
64 "using cubic interpolation instead"]);
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20567
diff changeset
65 der(:,1) = feval (func, z(1), u(:,1), args);
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20567
diff changeset
66 der(:,2) = feval (func, z(2), u(:,2), args);
20565
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
67 u_interp = hermite_cubic_interpolation (z, u, der, t);
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20567
diff changeset
68
20565
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
69 endswitch
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
70
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
71 endfunction
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
72
20847
ddc18b909ec7 codesprint: ode: first and second order interpolator for dense_output
jcorno <jacopo.corno@gmail.com>
parents: 20634
diff changeset
73 ## The function below can be used in an ODE solver to interpolate the solution
ddc18b909ec7 codesprint: ode: first and second order interpolator for dense_output
jcorno <jacopo.corno@gmail.com>
parents: 20634
diff changeset
74 ## at the time t_out using 2th order hermite interpolation.
ddc18b909ec7 codesprint: ode: first and second order interpolator for dense_output
jcorno <jacopo.corno@gmail.com>
parents: 20634
diff changeset
75 function x_out = quadratic_interpolation (t, x, der, t_out)
ddc18b909ec7 codesprint: ode: first and second order interpolator for dense_output
jcorno <jacopo.corno@gmail.com>
parents: 20634
diff changeset
76
ddc18b909ec7 codesprint: ode: first and second order interpolator for dense_output
jcorno <jacopo.corno@gmail.com>
parents: 20634
diff changeset
77 # coefficients of the parabola
ddc18b909ec7 codesprint: ode: first and second order interpolator for dense_output
jcorno <jacopo.corno@gmail.com>
parents: 20634
diff changeset
78 a = -(x(:,1) - x(:,2) - der(:).*(t(1)-t(2))) / (t(1) - t(2))^2;
ddc18b909ec7 codesprint: ode: first and second order interpolator for dense_output
jcorno <jacopo.corno@gmail.com>
parents: 20634
diff changeset
79 b = der(:) - 2*t(1).*a;
ddc18b909ec7 codesprint: ode: first and second order interpolator for dense_output
jcorno <jacopo.corno@gmail.com>
parents: 20634
diff changeset
80 c = x(:,1) - a*t(1)^2 - b*t(1);
ddc18b909ec7 codesprint: ode: first and second order interpolator for dense_output
jcorno <jacopo.corno@gmail.com>
parents: 20634
diff changeset
81
ddc18b909ec7 codesprint: ode: first and second order interpolator for dense_output
jcorno <jacopo.corno@gmail.com>
parents: 20634
diff changeset
82 # evauate in t_out
ddc18b909ec7 codesprint: ode: first and second order interpolator for dense_output
jcorno <jacopo.corno@gmail.com>
parents: 20634
diff changeset
83 x_out = a*t_out.^2 + b*t_out + c;
ddc18b909ec7 codesprint: ode: first and second order interpolator for dense_output
jcorno <jacopo.corno@gmail.com>
parents: 20634
diff changeset
84
ddc18b909ec7 codesprint: ode: first and second order interpolator for dense_output
jcorno <jacopo.corno@gmail.com>
parents: 20634
diff changeset
85 endfunction
20565
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
86
20903
3d3da166dac5 2015 Code Sprint: finish import of ode23 into core
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20847
diff changeset
87 ## The function below can be used in an ODE solver to interpolate the
3d3da166dac5 2015 Code Sprint: finish import of ode23 into core
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20847
diff changeset
88 ## solution at the time t_out using 4th order hermite interpolation.
20565
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
89 function x_out = hermite_quartic_interpolation (t, x, der, t_out)
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
90
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
91 persistent coefs_u_half = ...
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
92 [(6025192743/30085553152), 0, (51252292925/65400821598), ...
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
93 (-2691868925/45128329728), (187940372067/1594534317056), ...
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
94 (-1776094331/19743644256), (11237099/235043384)].';
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
95
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20567
diff changeset
96 ## 4th order approximation of y in t+dt/2 as proposed by Shampine in
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20567
diff changeset
97 ## Lawrence, Shampine, "Some Practical Runge-Kutta Formulas", 1986.
20565
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
98 dt = t(2) - t(1);
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
99 u_half = x(:,1) + (1/2) * dt * (der(:,1:7) * coefs_u_half);
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
100
20565
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
101 ## Rescale time on [0,1]
20567
ea6a1c00763a fix interpolation bug introduced with 87b557ee8e5d
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
102 s = (t_out - t(1)) / dt;
20565
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
103
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
104 ## Hermite basis functions
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
105 ## H0 = 1 - 11*s.^2 + 18*s.^3 - 8*s.^4;
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
106 ## H1 = s - 4*s.^2 + 5*s.^3 - 2*s.^4;
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
107 ## H2 = 16*s.^2 - 32*s.^3 + 16*s.^4;
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
108 ## H3 = - 5*s.^2 + 14*s.^3 - 8*s.^4;
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
109 ## H4 = s.^2 - 3*s.^3 + 2*s.^4;
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
110
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
111 x_out = zeros (rows (x), length (t_out));
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20567
diff changeset
112 x_out = (1 - 11*s.^2 + 18*s.^3 - 8*s.^4) .* x(:,1) + ...
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20567
diff changeset
113 ( s - 4*s.^2 + 5*s.^3 - 2*s.^4) .* (dt * der(:,1)) + ...
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20567
diff changeset
114 ( 16*s.^2 - 32*s.^3 + 16*s.^4) .* u_half + ...
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20567
diff changeset
115 ( - 5*s.^2 + 14*s.^3 - 8*s.^4) .* x(:,2) + ...
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20567
diff changeset
116 ( s.^2 - 3*s.^3 + 2*s.^4) .* (dt * der(:,end));
20565
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
117
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
118 endfunction
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20567
diff changeset
119
20903
3d3da166dac5 2015 Code Sprint: finish import of ode23 into core
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20847
diff changeset
120 ## The function below can be used in an ODE solver to interpolate the
3d3da166dac5 2015 Code Sprint: finish import of ode23 into core
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20847
diff changeset
121 ## solution at the time t_out using 3rd order hermite interpolation.
3d3da166dac5 2015 Code Sprint: finish import of ode23 into core
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20847
diff changeset
122 function x_out = hermite_cubic_interpolation (t, x, der, t_out)
3d3da166dac5 2015 Code Sprint: finish import of ode23 into core
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20847
diff changeset
123
3d3da166dac5 2015 Code Sprint: finish import of ode23 into core
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20847
diff changeset
124 s = (t_out - t(1)) / (t(2) - t(1));
3d3da166dac5 2015 Code Sprint: finish import of ode23 into core
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20847
diff changeset
125 x_out = zeros (size (x, 1), length (t_out));
3d3da166dac5 2015 Code Sprint: finish import of ode23 into core
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20847
diff changeset
126
3d3da166dac5 2015 Code Sprint: finish import of ode23 into core
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20847
diff changeset
127 for ii = 1:size (x, 1)
3d3da166dac5 2015 Code Sprint: finish import of ode23 into core
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20847
diff changeset
128 x_out(ii,:) = (1+2*s).*(1-s).^2*x(ii,1) + s.*(1-s).^2*(t(2)-t(1))*der(ii,1) ...
3d3da166dac5 2015 Code Sprint: finish import of ode23 into core
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20847
diff changeset
129 + (3-2*s).*s.^2*x(ii,2) + (s-1).*s.^2*(t(2)-t(1))*der(ii,2);
3d3da166dac5 2015 Code Sprint: finish import of ode23 into core
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20847
diff changeset
130 endfor
3d3da166dac5 2015 Code Sprint: finish import of ode23 into core
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20847
diff changeset
131
3d3da166dac5 2015 Code Sprint: finish import of ode23 into core
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20847
diff changeset
132 endfunction