Search All of the Math Forum:
Views expressed in these public forums are not endorsed by
Drexel University 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 pre-allocated 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 re-calculate 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 pre-allocating 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, pre-allocate 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 pre-allocate 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
|
|
|
|
