Search All of the Math Forum:
Views expressed in these public forums are not endorsed by
NCTM or The Math Forum.



Re: ode45 solver and outputfcn
Posted:
Apr 26, 2011 2:12 PM


On 26 Apr, 00:08, "Nasser M. Abbasi" <n...@12000.org> wrote: > On 4/25/2011 2:19 PM, Allamarein wrote: > > > > > > > > > > > I often solved my ODEs by the sintax: > > > [T,Y] = ode15s(@vdp1000,[0 3000],[2 0]); > > > and let's say the problem was: > > > function dy = vdp1000(t,y) > > dy = zeros(2,1); % a column vector > > dy(1) = y(2); > > dy(2) = 1000*(1  y(1)^2)*y(2)  y(1); > > > it works fine. > > > Now let's suppose my ODEs is > > > function dy = vdp1000(t,y) > > dy = zeros(2,1); % a column vector > > % K the new value > > K = newvalue(t,y) > > dy(1) = y(2); > > dy(2) = 1000*(1  y(1)^2)*y(2)  y(1) + K; > > > It means there is a K value computed at each time by a function called > > newvalue. > > I would have this value T in output at each output time. For example: > > [T,Y, K] = ode15s(@vdp1000,[0 3000],[2 0]); > > > Untill now, at the end of my simulation, I "rebuild" T, but it is a > > time consuming procedure and dumb. > > I mean that: > > [T,Y] = ode15s(@vdp1000,[0 3000],[2 0]); > > for i=1:length(T) > > K(i) = newvalue(T(i), Y(i,:) ) > > end > > > I suspect the option 'outputfcn ' could help me. > > Someone knows a better method? > > The problem is that we do not know how many time > steps the ode solver will use, else you could > have preallocated a vector large enough to hold the 'k' > values, and simply populated this vector as the program runs, > so at the end you do not need to recalculate it. > > These solvers are adaptive, and even if one gives it a > time vector as in > > timeSpan = [t0:delt:tFinal] > > This does _not_ mean that it will make only one time > step each delt, as it can make more within that delt > time. delt acts as a worst case limit. > > The ode solvers decide on the step size, but delt is taken > as an upper limit. This can be easily seen: > >  > function foo > > delT = 0.1; > timeSpan = 0:delT:1; > counter = 0; > > [t,y]=ode45(@RHS,timeSpan,0); > > fprintf('length of time vector = %d, number of time ODE called %d\n',... > length(t), counter); > > function dydt=RHS(t,y) > counter = counter + 1; > dydt=exp(t); > end > end >  > > The above program generates this output: > > length of time vector = 11, number of time ODE called 85 > > So you see, even if the time step was 0.1, ode made 85 > steps. Much more than we expected. > > So preallocating vector for K will not help. > Second option is to simply append and entry to the > k vector. Will be slower ofcourse, but you would > not have to recalculate everything again. You need > to decide which is more efficient: > >  > function foo > > delT = 0.1; > timeSpan = 0:delT:1; > k = zeros(length(timeSpan)); > counter = 0; > > [t,y]=ode45(@RHS,timeSpan,0); > > fprintf('length of time vector = %d, number of time ODE called %d\n',... > length(t), counter); > > k; %now k has all the calculated values > > function dydt=RHS(t,y) > counter = counter + 1; > k(counter) = t^2; %or whatever. WARNING Dynamic allocation possible > dydt=exp(t); > end > end >  > > I am not sure what @output will help with here. > > A third option is this: Run the ode solver you have, on the > same problem as a 'dry' run first, find the number of steps > used by the solver. Next, preallocate a k vector of this size, > and now you can avoid the dynamic allocation issue. Make > sure you change nothing in the problem parameters so > that the number of steps actually used by ODE do not change. > > A fourth option, is to write your own solver (say RK4) > directly (easy to do) and this way you have control on > the number of step made, and can preallocate things to > the exact sizes. But you would lose the adaptive step > size part. This can be important, for say stiff problems > and such. > > Nasser
Thanks Nasser. Your suggestions allow to get around the problem. I would solve my problem interfacing with odepackage. In my opinion this is the most correct method and it would permit me to save a lot of time



