ODE45 errer :input argument y is undefined

hi everybody, i've got a problem and it drives me crazy...

In a function i defined a systeme :

function  dy=f_iso(t,y)

global     Q b K n E sigma0 epoint

if abs(y(1))-(Q*(1-exp(-b*y(3))))-sigma0<0

      dy(3)=0;

else

      dy(3)=((abs(y(1))-(Q*(1-exp(-b*y(3))))-sigma0)/K)^n;

end

dy(2)=sign(y(1))*dy(3);

dy(1)=E*(epoint-dy(2));

dy=dy';

end 

Then i've used ode45 such as : global  Q b  n E sigma0 epoint

Q=-150;

b=5;

n=6;

sigma0=200;

E=140000;

epoint=0.001;

[t,y]=ode45('f_iso',[0 10],[0 0 0]);

sigma1=y(1,:);

ep1=y(3,:);

e1=sigma1/E+ep1;

plot(e1,sigma1)

but there is qn error :

Input argument "y" is undefined.

Error in ==> f_iso at 8

and i don't know whyt...

thanks

Answer 1

[Jan]

Please post a copy of the complete error message and mention the line, which causes the error. If "line 8" is dy(2)=sign(y(1))*dy(3);, this would be very strange, because the variable y is used in if abs(y(1))-(Q*(1-exp(-b*y(3))))-sigma0<0 already.

Therefore I guess, some standard bugs are responsible. Did you shadow a built-in function by a script, which contains a clear all? To check this, set a breakpoint in the function f_iso and step through the code line by line trying to jump into subfunctions, e.g. by the F11 key.

Btw., the function to be integrated must be smooth, otherwise the step-size control of ODE45 can lead to unexpected effects:

1. ODE45 integrates right over the discontinuity without noticing this. Then the result has a poor accuracy only.
2. ODE45 finds the discontinuity and stops with a warning message "Unable to meet integration tolerances without reducing the step size below the smallest value allowed".
3. ODE45 reduces the step size to such a tiny value, that the integration takes hours to run and the accumulated rounding error dominates the solution. This is exactly what the step size control should avoid.

Which of these cases happens can depend on minimal variations of the initial values or parameters, and therefore this must be treated as a bug. Use an event function to toggle values of a parameter.

[EDITED] Some references:

• http://saba.kntu.ac.ir/eecd/ecourses/Matlab%20premier/Chapter08.pdf
• http://comjnl.oxfordjournals.org/content/13/4/401.full.pdf
• http://pcmap.unizar.es/numerico/upload/DisconRep-2002.pdf

The last one explains a method to let a DOPRI45 detect and handle a discontinuity. It would be a nice enhancement for Matlab's ODE45. It could be expanded to handle parameter changes triggered by event functions without the need to restart the integration.

[Srikumar Gopalakrishnan]

Hello Jan Can you please check the third link it would be of great help if you could upload the document in the third link. There is an error in this link (<http://pcmap.unizar.es/numerico/upload/DisconRep-2002.pdf>)

[D Zoff]

Hi Jan, the third link doesn't seem working now. Do you know where I can find this pdf document?

Also, may I please ask, if I suspect that my function to be integrated contains discontinuities, how can I detect it and confirm that it does contain discontinuities. My ode45 runs very slow at a certain time point so I think my function might have discontinuities. I am studying the flying trajectories of atoms in a molecule due to strong interatomic forces. My ode function is to calculate interatomic forces (not obtained via Coulomb's law but from results of electronic structure calculation, so the forces are not a function of any obvious quantities, which I guess it is where the discontinuities from) on every atoms in a molecule at a time point, and so acceleration on each atom is then calculated by F=ma, and the velocity can be obtained, and the new location of atoms can be found. My initial condition is zero velocity on all atoms, a given initial locations on all atoms and very strong interatomic forces on all atoms. I integrated this ode function over a certain time length (about 100 femtosecond). My problem now is that, I get integration stuck somewhere around 8 femtosecond and the time increment is really small. Is there anything I need to change?

[Jan]

I do not know where this paper can be found today. Perhaps this helps: https://www.researchgate.net/publication/223157481_Efficient_integration_over_discontinuities_in_ordinary_differential_equation_simulation

You can simply see in the code of the function to be integrated, if it contains discontinuities: Look for if, max, min, abs, tan (with arguments at pi*(2N-1)/2), linear interpolations, rand, etc.

A small time increment can have other reasons, e.g. high frequency changes of the forces. The problem might be stiff also.

[D Zoff]

Thank you so much, Jan. I will read through the paper you mentioned.

In fact, what I experienced is, at first 10 femtosecond integration time, it all looks fine. As soon as t passed 10 femtosecond, I noticed the time went back to 8.5 femtosecond and integrated again but very slowly from 7.5 ft to 8.2 ft and almost stayed at 8.2 ft for hours. What could the possible reason for this went-back-and-restart integration behaviour?

I did try to take the results when t=10 ft (last timepoint before t went back) and use it as a new initial condition for the rest integration time (i.e. 10fs-100fs), the ode then successfully terminated at 100fs. It seems to me that if I can somehow manage to 'force' my integration time t not go back and restart, I can get what I want. But it is a reasonable strategy? If so, how can I do it?

[Jan]

Your function might have a pole. Then a step is tried, but the difference between the two embedded integration schemes of ODE45 is too large to be accepted. We could examine your function, if you post it.

[D Zoff]

Here is the simplified version of my ode funtion.

The quantum chemistry software is an external program that I can call and manipulate its result log file from within matlab. I only need it to give me interatomic forces, which is used to determine accelerations.

I used this ode function to look at molecules at different charge states. Presumably, the higher the charge state, the larger the interatomic repulsive forces. For many different charge states, it works fine (speedy finish time and reasonable results). Only for a certain charge state, it has this issue.

Thank you very much.

My program:

tspan=[0 1e-13];

[t,y] = ode45(@NEfun,tspan,yinit);

function dydt = NEfun(t,y)

for k=1:3*numberofatoms

      dydt(k)=y(3*numberofatoms+k);   % pass velocities from y's last numberofatoms*3 elements to dydt's first numberofatoms*3 elements

end

geo=zeros(numberofatoms, 3);

for k=1:numberofatoms

    geo(k,1)=y(1+3*(k-1));     % pass atomic coordinates from y's first numberofatoms*3 elements to matrix geo
    geo(k,2)=y(2+3*(k-1));
    geo(k,3)=y(3+3*(k-1));

end;

%then:

%1) construct a quantum chemistry software input file by using the matrix geo as the molecular geometry

%2) execute the quantum chemistry calculation

%3) when the quantum chemistry calculation finished, a log file containing interatomic forces is produced.

%4) extract the interatomic forces from the log file and store it into a matlab matrix called force (size: numberofatoms X 3)

for n = 1:numberofatoms

      for l=1:3

          dydt(3*numberofatoms+3*(n-1)+l) = force(n,l)/Mass(n);  %calculate and pass accelerations to dydt's last numberofatoms*3 elements

      end

end

dydt=dydt';  

%=========list  of variables=========

t             % evaluation points

y             % a numberofatoms*3*2 long solution vector that first numberofatoms*3
              % elements are coordinates of atoms and the last numberofatoms*3
              % elements are velocities of atoms.

dydt          % a numberofatoms*3*2 long vector that first numberofatoms*3 
              % elements are velocities of atoms and the last 
              % numberofatoms*3 elements are accelerations of atoms

tspan         % interval of integration time, set to from zero to 100 fs

numberofatoms % Number of atoms, positive integer

geo           % a numberofatoms X 3 matrix containing x y z coordinates of all atoms

force         % a numberofatoms X 3 matrix containing forces of all atoms (in x y z directions)

mass          % a numberofatoms long vector containing masses of all atoms

%====================================

[D Zoff]

So the basic idea of this ode is, at each time step, use current atomic coordinates to perform a quantum chemistry calculation to get interatomic forces, then to get accerleration from F=m*a, then since current coordinates and velocities are known, ode can tell me the next step coordinates, which can then be used for a new quantum chemistry calculation to get a new force at next time step, and so on.

[D Zoff]

Hi Jan,

Sorry if the program I wrote above is too messy. I am going to make it simpler if that helps.

1) Problem 1

If say, I am trying to find the trajectories evolution (and velocity evolution) of two objects (say atoms) that possess electrical charges q1 and q2, which makes them repulsively flying apart due to the Coulomb’s law. So the initial condition is that: they are resting at two different locations x10, y10, z10 and x20, y20, z20, and they both have zero velocity. In terms of programming, y0= [x10 y10 z10 x20 y20 z20 0 0 0 0 0 0 ] (i.e. coordinates of the initial location of two objects followed by coordinates of the initial velocities of the two objects).

Inside the function to be solved, dydt=myfunction(t,y), I pass the initial velocities (i.e. the last 6 elements in y0) to the first 6 elements of dydt. i.e. dydt(1st 6 elements) = y(last 6 elements). Then I calculated the forces acted on the two objects by using Coulomb’s Law, F=k x q1q2/r^2. With known mass m, I can use Newton’s equation, F= m x a to calculate the acceleration a, and I pass the acceleration to the last 6 elements of dydt. So the dydt will be something like = [velx1 vely1 velz1 velx2 vely2 velz2 accx1 accy1 accz1 accx2 accy2 accz2]

For this problem, since the Coulomb Forces will be smoothly decreasing (so do the velocities), it should be handled perfectly by Matlab’s ode45 without theoretical flaw, am I right?

2) Problem 2

My real problem is similar to problem 1 stated above, but with an annoying modification.

Inside the function to be solved, I do not use Coulomb’s law to calculate the interatomic forces anymore. Instead, I entirely rely on an external software to tell me what the interatomic forces will be, at every single time step. That means, inside the function to be solved, I wrote some lines to send out the coordinates of the two atoms to an external software, and ask that software to send me back the respective interatomic forces. These interatomic forces do depend on the coordinates of the two objects, but they also depend on many other factors which are outside the scope of the variables we know here (mostly related to electronic structures of atoms which are handled by that software). These new interatomic forces are close to the forces calculated by using Coulomb’s Law in Problem 1. To some extend, it can be thought of as a value that is slightly deviated from the interatomic forces calculated by Coulomb’s Law in Problem 1, but the magnitude of deviation is unknown, and is not a function of anything variables we know here.

I have encountered many issues when trying to run problem 2:

1) Is the function described in problem 2 a non-smooth function with discontinuity?

2) If it is, how should I handle it? I don’t quite see that there is an event that I can monitor to stop and restart the integration process. To my understanding, every time I got a force value from the external software, it introduced a discontinuity, am I right?

3) The range of the solution y is hugely different between coordinates (1st 6 elements) and velocities (last 6 elements). Both atoms have coordinates in the range of about 0 ~ +/- 100, but have velocities in the range of about 0 ~ +/- 10000. I guess the error tolerance should be specified differently. So I gave something like [1e-6 1e-6 1e-6 100 100 100] to Abs Tol, and I kept the default Rel Tol 1e-3 unchanged. Is this sensible? I can only think of the fact that Abs Tol for coordinates should be much smaller than Abs Tol for velocities, but how to determine the most suitable values? I did noticed that when I first ran problem 2 with default Abs Tol, i.e. 1e-6 for both coordinates and velocities, my ode was stuck at some time step and moving extremely slowly.

4) I added a few fprintf commands inside the function to be solved to write the latest data (time step, coordinates, velocities, forces) to a text file at each time step, so that even my program crushes at any point, I still have the data I wanted. When finished one whole ODE. I noticed that the time steps in the given output t is far fewer than the actual time steps I recorded from fprintf commands throughout the whole process (is this because the actual time steps are the ‘tried’ time step?). Of course, time steps in t do exist in the actual time steps. Also, the actual time steps repeat itself once at every 6 time steps. Sometime, after the repeated time steps, it went back to an earlier time step to try again! Why there are two different time steps? Why the actual time steps decided to go back at some point? Because the process of asking the external software to get me a force value is the most time consuming part of the each circle (about 10 min each time. So if there is a 1000 actual time steps, it will be about 7 days needed!), what is the correct strategy to significantly shorten the total calculation time, while maintaining acceptable errors? The time that the external software needs for each calculation cannot be shorten, so the only way is to have fewer steps.

5) After using ode45, I did try using ode15. But it was even worse as the integration spent a huge number of time steps at 0, and only jump to tfinal at the very last few time steps. Besides, the total number of time steps is not reduced, meaning that there is not much benefit.

I was trying the make a ‘simpler’ description, but somehow inevitably ended up with such a long paragraph. I do apologise for that.

Answer 2

[tony]

my program me finally work, but i don't why...

[Alioua Abderahmane]

hello!

i would like to ask you that i have problem with ode45 in line 195 knowing that i used the same things in the prevididn't receive any error

Undefined function or variable 'Prav_Nabl_NL_2mayat_1a'.

Error in Two_Inv_Pend_Nabl_Reg_070620_al>@(t,x)Prav_Nabl_NL_2mayat_1a(t,x,K1,L1,C)

Error in odearguments (line 87)

f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)

odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in Two_Inv_Pend_Nabl_Reg_070620_al (line 196)

[t,x] = ode45(@(t,x) Prav_Nabl_NL_2mayat_1a(t,x,K1,L1,C),[t0 tk],x0);