MATLAB Answers

ODE45 errer :input argument y is undefined

128 views (last 30 days)
tony
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

Answers (2)

Jan
Jan on 22 Jan 2013
Edited: 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:
  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:
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.
  8 Comments
D Zoff
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.


tony
tony on 22 Jan 2013
my program me finally work, but i don't why...
  1 Comment
Alioua Abderahmane
Alioua Abderahmane on 10 Jun 2020
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);

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!