How can I save one of variable in my function ?

Hi,

I want to save one of the time-dependent variables ( "NOD1") given in FUNCTION.m , yet I could not figure it out. Can You Please help me with that?

I have attached here the FUNCTION.m and MAIN.m

I want to save all results of "NOD1".

Is that possible?

Thanks for your help.

Kind regards,

[James Tursa]

Looks like your derivative function is not continuous (i.e., has an if-then-else part that causes the output to jump at certain points). This may cause problems with ode45( ).

Answer 1

[Torsten]

After ode45 has finished, just recalculate NOD1 in the calling program for the time vector t:

[t,x]=ode45(@odeHar,tspan,x0);

C_n = 300*10^3;
NOCo1 = x(:,1).*x(:,3);
NOD1 = C_n*x(:,3).*(NOCo1 < 0) + 100*x(:,3).*(NOCo1 >= 0)

[Walter Roberson]

This approach makes a lot more sense. Calculate what the NOD1 would be for the actual points that are output, giving you one NOD1 for each actual output, in a known order. These would be meaningful NOD1 values. Recording each NOD1 value as odeHar is called would not give you meaningful NOD1 values.

[ADNAN]

many thanks @Walter Roberson and @Torsten

Answer 2

[Walter Roberson]

Don't do that.

ode45 works by evaluating the given function 6 times per trial, with the inputs carefully chosen according to the current boundary conditions and a current internal step size. One of the evaluations is a cross-check on the others. If the cross-check fails then ode45 re-uses the same initial position with a smaller step size. A number of cross-checks in a row can fail... I estimate 40 or more in a row could fail before it gave up in some cases.

You are asking that all those NOD1 be recorded, even though many in a row might not be used in the final position.

If you were walking in the dark on hills with a wind-up strobe flashlight, would you be wanting to record all of the heights or only the heights corresponding to the direction you ended up travelling?

Now, in the case where the cross-check succeed, then with the default options, ode45 not only records one of the outputs: it also interpolates three intermediate positions between the previous position and the current position, without calling the function. So there are outputs for which no corresponding NOD1 was ever calculated. So not only would there be many useless NOD1 values recorded, there would be outputs for which no NOD1 was recorded!

Recording "all" values of an intermediate for ode45 is almost always a bad idea.

[Walter Roberson]

 z_t=interp1(tspan,acc,t);

ode45 is not a suitable tool for cases in which there are non-smooth forces being applied, such as external acceleration, or such as regular injections of drugs, or such as bouncing off of something. The mathematics of the ode45 computation is invalid in such cases. You must not cross a non-smooth boundary within any one call to ode45.

In the special case that the external force can be understood to be a sampling of a continuous but variable function, then using interp1 with spline interpolation can be mathematically valid.

[ADNAN]

Thanks for that, then what is the solution? I need to save "NOD1". thanks again.

[Walter Roberson]

if NOCo1<0

You need a terminal event function to handle that case: you need to tell ode45 to stop integrating at the boundary where NOCo1 becomes 0, and then call ode45 again to continue from where it left off.

[Walter Roberson]

I need to save "NOD1"

Suppose you have the best case, where the step succeeds. 6 function calls were made, and you do not record the associated parameters and you would have to dig into the code to find out how the parameters are chosen, so you would have 6 recordings of NOD1 with at most one of them having any relationship to a recorded set of boundary conditions. And then it would have interpolated three output boundary conditions without having called the function. So in the best case, you would have repeated stretches of

[interpolated boundary with no NOD1 calculated, interpolated boundary with no NOD1 calculated, interpolated boundary with no NOD1 calculated, non-interpolated boundary with 6 NOD1 calculated]

with it being unspecified which of the 6 NOD1 values would correspond to the non-interpolated boundary.

I said "by default" but I just noticed you are not using the defaults that applies to: you are using a vector of more than two time steps. When you use a vector of time steps, all output boundary conditions (except possibly the first and the last) are interpolated when ode45 detects that it has crossed the requested time, and NOD1 would not be recorded for any of those outputs.

You might as well use a random number generated for all of the good it is going to do you to have those NOD1 values recorded.

[ADNAN]

thanks for that @Walter Roberson