Asked by tony
on 22 Jan 2013

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 by Jan
on 22 Jan 2013

Edited by Jan
on 16 Apr 2015

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:

- ODE45 integrates right over the discontinuity without noticing this. Then the result has a poor accuracy only.
- 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".
- 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.

D Zoff
on 5 May 2017

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
on 5 May 2017

D Zoff
on 31 Aug 2017

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.

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.