@pagan A bit inconsistent.

The code does not work for getting a function of the solution like F[i]=F[i-1]+f(sol[i])...

@pagan Well, I think this is what I meant.

However, I meant not F[i]=sol(y), but F[i]=f(sol(y)).

Say, F[i]=F[i-1]+sol[i](y).

I hope this will not create any difference.

And yes, without rmemebr option this is rather slow.

I am afraid this will kill the idea, since in reality I need to implement this to a system of 6 DE's...

we will see.

Anyway, this is the closest thing I've got since the beginning  of asking questions here.

I guess this should be like

restart:

sys:= diff(y(t),t) = cos( 1/2*y(t)-x(t)-Pi^2 ):
ini:= y(0)=10*sqrt(2)/2:

F[0]:=sin:

N:=10:

for i from 1 to N do
   newsys:=subs(x=F[i-1],sys);
   #print(newsys);
   sol[i]:=dsolve([newsys, y(0)=10*sqrt(2)/2], numeric,
                  output=listprocedure, known=F[i-1]):
   eval(y(t),sol[i])(0.2);
   #F[i]:=unapply('evalf'(subs(y(t)=eval(y(t),sol[i])(t),x=F[i-1],rhs(sys))),
   #              t,numeric,proc_options=[remember]):
   F[i]:=F[i-1]+subsop(3=remember,eval(y(t),sol[i]));
end do:

@pagan Well, I think this is what I meant.

However, I meant not F[i]=sol(y), but F[i]=f(sol(y)).

Say, F[i]=F[i-1]+sol[i](y).

I hope this will not create any difference.

And yes, without rmemebr option this is rather slow.

I am afraid this will kill the idea, since in reality I need to implement this to a system of 6 DE's...

we will see.

Anyway, this is the closest thing I've got since the beginning  of asking questions here.

I guess this should be like

restart:

sys:= diff(y(t),t) = cos( 1/2*y(t)-x(t)-Pi^2 ):
ini:= y(0)=10*sqrt(2)/2:

F[0]:=sin:

N:=10:

for i from 1 to N do
   newsys:=subs(x=F[i-1],sys);
   #print(newsys);
   sol[i]:=dsolve([newsys, y(0)=10*sqrt(2)/2], numeric,
                  output=listprocedure, known=F[i-1]):
   eval(y(t),sol[i])(0.2);
   #F[i]:=unapply('evalf'(subs(y(t)=eval(y(t),sol[i])(t),x=F[i-1],rhs(sys))),
   #              t,numeric,proc_options=[remember]):
   F[i]:=F[i-1]+subsop(3=remember,eval(y(t),sol[i]));
end do:

Thanks for help.

It seems this is what I need.

However the line

F[i]:=unapply('evalf'(subs(y(t)=eval(y(t),sol[i])(t),x=F[i-1],rhs(sys))),
                 t,numeric,proc_options=[remember]):

means that F[i] You substitute the calculated solution into the very same DE, am I right?
Is it possible to do the same with some other function and not with the system sys itself?
Say, i need F[i] to be the function of solution to sys itself.
Again thanks

Thanks for help.

It seems this is what I need.

However the line

F[i]:=unapply('evalf'(subs(y(t)=eval(y(t),sol[i])(t),x=F[i-1],rhs(sys))),
                 t,numeric,proc_options=[remember]):

means that F[i] You substitute the calculated solution into the very same DE, am I right?
Is it possible to do the same with some other function and not with the system sys itself?
Say, i need F[i] to be the function of solution to sys itself.
Again thanks

@pagan Right;

However both variants do not solve my problem, since I need the analytical form of the function explictly as function of time.

To this end the piecewise option is better, but since I use this function in further computatios this does not help either.

I think I should proceed with polynomila interpolation of a solution structure given as an Array.

Hopefully this will give me some polynomial and not piecewise function, which may be then will be simpler...

@pagan Right;

However both variants do not solve my problem, since I need the analytical form of the function explictly as function of time.

To this end the piecewise option is better, but since I use this function in further computatios this does not help either.

I think I should proceed with polynomila interpolation of a solution structure given as an Array.

Hopefully this will give me some polynomial and not piecewise function, which may be then will be simpler...

Thanks for the answer.

Basically I would like to know how to extract those data from the matrix F[2,1].

How one can do that?

Sorry for being not precise

Thanks for the answer.

Basically I would like to know how to extract those data from the matrix F[2,1].

How one can do that?

Sorry for being not precise

Thanks a lot, this is exactly what I have been looking for

Thanks a lot, this is exactly what I have been looking for

With the progress so far I may reformulate the question more specifically:

How one may evaluate functions of solutions of DE system to a function of time only?

With one variable the composition operator works, namely

x(y,t) with y=proc(t)..proc may be evaluated to x(t) by using (x@y)(t) if I understood correctly.

However, for function of several variables it turns to be harder...

Say, given u(x,y,t) and solutions x=proc(t)..t, y=proc(t)..t may I define u(t)?

Thanks,

I will try it.

However I must admit I do not fully yet understand the usage of this option.

Namely, I would like to evaluate a function u(x(t),y(t),t) as a funxtion u(t) for given solution x(t), y(t).

I am not sure the remember option would help here...

Thanks,

I will try it.

However I must admit I do not fully yet understand the usage of this option.

Namely, I would like to evaluate a function u(x(t),y(t),t) as a funxtion u(t) for given solution x(t), y(t).

I am not sure the remember option would help here...

Thanks for the answer.

I speak Russian but do not have a Russian keyboard with me right now.

However: The problem is NOT in storing the results, but in the way optimal control is defined.

I have control u=f(x,y,t) with x and y state and costate variables respectively. I need to evaluate u as f(t) at every step to avoid computational burden.

For example  I have u(x,y,t)=x^2+y^2-ln(t) and x=f(t), y=f(t), but defined through numerical procedures as proc(t)... proc.

I need maple to evaluate u(t) given this. However with procedure-defined functions this seems to be difficult...

This control function is supposed to be used then for the next step of the procedure, namely

0. Take initials;

1. Compute x,y as functions from initials and t

2. Use x, y to define u(t) - control function

3. Update x,y from 1.

And I have the problem at step 2.