bugGNU Octave - Bugs: bug #32959, lsode use "1/f" instead...

 
 

bug #32959: lsode use "1/f" instead "f"

Submitter:  None
Submitted:  Wed 30 Mar 2011 09:28:49 PM UTC
   
 
Category:  Libraries Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Works For Me Assigned to:  None
Originator Name:  Valber Originator Email:  -email is unavailable-
Open/Closed:  * Closed Release:  * 3.3.50
Operating System:  * GNU/Linux Fixed Release:  None
Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Sun 10 Apr 2011 02:19:31 PM UTC, comment #7: 

There is no mistake here.  "bdf" and "stiff" are the same solver.  Likewise "adams" and "non-stiff" are the same solver.  We just give people the option of using either name.

Rik <rik5>
Group administrator
Sat 09 Apr 2011 07:12:58 PM UTC, comment #6: 

I sorry, and thank you for help. I understand you.
I translate Octave doc into russian.

Another mistake, I open the report??
listing (Octave 3.2.4)

>>> lsode_options("integration method")

ans = stiff

>>> lsode_options("integration method","bdf")
>>> lsode_options("integration method")

ans = stiff

>>> lsode_options("integration method","adams")
>>> lsode_options("integration method")

ans = non-stiff

>>> lsode_options("integration method","stiff")
>>> lsode_options("integration method")

ans = stiff

Valber <psi13art>
Fri 08 Apr 2011 10:47:18 PM UTC, comment #5: 

The problem is the order of your arguments in the function which calculates the derivative.  According to the lsode documentation ('help lsode'):


The first argument, FCN, is a string, inline, or function handle
that names the function f to call to compute the vector of right
hand sides for the set of equations.  The function must have the
form

     XDOT = f (X, T)

in which XDOT and X are vectors and T is a scalar.


In your sample problem you are using Y instead of X as the dependent variable and X instead of T as the independent variable.  Therefore, the functional form should be

YDOT = f (Y, X)


However, without specific instructions the inline() function orders its inputs alphabetically.  For example,

f = inline ("-x/y")
f = f(x, y) = -x/y
      ^ NOTE that x is first, and you require f (y, x) according to lsode

# Specify order of arguments
f = inline ("-x/y", "y", "x")
f = f(y, x) = -x/y
      ^ NOTE y is now first.

# This code now works using 'f' above
x = -5:0.1:5;
y = lsode (f, 5, x);
plot (x,y);


Of course, you can still create a problem with


f = inline ("-y/x", "y", "x");


But, this shouldn't be a surprise because the function that calculates your derivative (y/x) blows up in the neighborhood of x = 0.  You can see this by using a range which doesn't go over the discontinuity.


f = inline ("-y/x", "y", "x");
x = -5:0.1:-0.1;
y = lsode (f,5,x);
plot (x,y);



Rik <rik5>
Group administrator
Fri 08 Apr 2011 01:56:52 PM UTC, comment #4: 

Ok, I'm try following example on 3.2.4

listing
=====================
GNU Octave, version 3.2.4
Copyright (C) 2009 John W. Eaton and others.

octave:1> x=-5:0.1:5;
octave:2> f=inline("y/x");
octave:3> lsode(f,5,x);
octave:4> plot(x,ans) %Ok, it is line!
octave:5> f=inline("-x/y");
octave:6> lsode(f,5,x);
 LSODE--  AT T (=R1) AND STEP SIZE H (=R2), THE   
       CORRECTOR CONVERGENCE FAILED REPEATEDLY    
       OR WITH ABS(H) = HMIN  
      In above,  R1 = -0.3729451955262-153   R2 =  0.1137126019255-173
error: lsode: repeated convergence failures (t = -3.72945e-154perhaps bad jacobian supplied or wrong choice of integration method or tolerances)

%Where is my Arc???
octave:6> f=inline("-y/x");
octave:7> lsode(f,5,x);
octave:8> plot(x,ans) %plot arc
=====================

Valber <psi13art>
Thu 07 Apr 2011 11:26:20 PM UTC, comment #3: 

The version you are using, 3.3.50, is a point release for developers.  It may have errors.

I ran your sample code without problems on the two stable released versions of Octave (3.2.4 and 3.4.0).

I would suggest changing to a stable version of Octave and re-opening this issue if you are still unable to get your code working.

Rik <rik5>
Group administrator
Thu 07 Apr 2011 09:08:59 PM UTC, comment #2: 


>value of the vector x

x=-5:0.1:5;

 solution equation  dy/dx=y/x => y=constant*x
 solution equation  dy/dx=-x/y => y=sqrt(const^2 - x^2) or y=-sqrt(const^2 - x^2)

If I type

x=-5:0.1:5;
f=inline("y/x");
lsode(f,5,x);

Octave make Error
But, if I type

x=-5:0.1:5;
f=inline("x/y");
lsode(f,5,x);
plot(x,ans) %plotting Line y=constant*x

=>LSODE using 1/f function instead f

Valber <psi13art>
Thu 31 Mar 2011 09:31:07 PM UTC, comment #1: 

There isn't enough information here to debug the problem.  What is the value of the vector x (the points where you want the function integrated)?

Also, your say that "y/x" is the problem, but your test results show "-x/y" being the problem.

Rik <rik5>
Group administrator
Wed 30 Mar 2011 09:28:49 PM UTC, original submission:  

Solve Ordinary Differential Equations
     dy             
     -- = f (y, x)
     dx
I set inline function
f=inline("y/x")

and lsode make eror

if set f=inline("x/y")

output is correct

TEST:

>>> f=inline("-x/y");
>>> lsode(f,5,x);

error: lsode: repeated convergence failures (t = -3.72945e-154perhaps bad jacob
ian supplied or wrong choice of integration method or tolerances)

>>> f=inline("x/y");
>>> lsode(f,3,x);
>>> plot(x,ans)
>>> f=inline("-y/x")
>>> lsode(f,3,x);
>>> plot(x,ans)

Anonymous

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

No files currently attached

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by psi13art (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by None (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 3 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2011-04-07 rik5 StatusNeed Info Works For Me
        Open/ClosedOpen Closed
    2011-03-31 rik5 StatusNone Need Info

    Back to the top

    Powered by Savane 3.13-f8d8.
    Corresponding source code