Structure of code for ODE solution in Octave

[Dr.D]

TL;DR

Is there a way to do preliminary calculations required only one time and then pass the results (via the parameter list?) to the derivative function?

The problem of interest at the moment is the solution of a simple damped oscillator problem,
xddot+2*zeta*wn*xdot+wn^2*x = 0
Everything seems to work fine provided I include the values of zeta and wn inside the function that defines the derivative values.

But zeta and wn are actually calculated from the mass, damping, and stiffness values, M, B, and K, as
wn=sqrt(K/M)
zeta=B/(2*wn*M)
What I'm seeing in Octave is that these calculations must be included inside the derivative evaluation function and executed repeatedly with every function call. Is that true?

For this simple problem this is a trivial waste, but in a more complex problem, with many parameters, it would be desirable to calculate items like wn and zeta before the first call to the function and then pass the computed values to the function as arguments. Is this possible? It does not seem to work, and every example I've found so far, does all the parameter evaluation within the derivative function.

Here is an example (that does not quite work!) of the question above:

% Test7.m
1;     % a do-nothing line to avoid file being a function definition
function zdot=g(z,t)
  M=3;                                 % mass
  K=9;                                  % stiffness
  B=.6;                                 % damping coeff
  wn=sqrt(K/M);                % undamped nat'l freq
  zeta=B/(2*M*wn);         % damping factor
  zdot(1)=z(2);                                              % x-dot
  zdot(2)=-2*zeta*wn*z(2)-wn^2*z(1);   % x-ddot
endfunction
tf=2;                                  % final time
soln=lsode("g",[1, 0],t=linspace(0,tf,100));
clear; clf; hold on;
x=soln(:,1)
xd=soln(:,2);
plot(t,x,'k',t,xd,'b');                % plot x(t) and xdot(t)
plot([0 tf],[0 0],'k')                 % draw time axis

In the context of the question above, is it not possible to move lines 4-8 outside the function definition?

As an additional question, does anyone see why this code will not run? If I remove the semi-colon at the end of the call to lsode, I get a two column output that appears to be soln(100,2), but this matrix never shows up in the Workspace list. I'm baffled!

PS: OK, I got the code to function correctly. The 'clear' after the lsode call was wiping out the solution. But the central question of avoiding inclusion of the parameter calcs remains.

 

[jedishrfu]

Does octave have a global keyword where inside the function you could say

Global M, K, B;

And pull those assignment statements outside the function g() ?

https://octave.org/doc/v4.2.1/Global-Variables.html

 

[Dr.D]

Thank you, @jedishrfu. Yes, there is a global declaration available in Octave, and I can make it work that way. But something seems to not work correctly even done that way.

If I use global variables, they have to be first declared outside the derivative routine, and then named again inside the derivative routine to make the connection. I've done this, and it works. But ...

If I physically modify the code byretype the line setting a global variable, it does not change if the code has been previously run. The old value is "stuck" in the machine. The only way I have found to this point to change a parameter value is to modify the code setting the value, save the code, exit Octave, restart Octave, reload the code, and then execute with the new values. Surely (I hope!) there has to be a better way!

 

[jedishrfu]

Odd it seems it should work unless you need to use global in both places. When you get this type of behavior it may mean a local variable of the same name has been defined and is hiding the global one.

 

[Dr.D]

@ jedishrfu: Thanks for your help.
From what I've found on the 'net, apparently Octave requires that global variables be tagged as such in the main program and again in each function using them. It seems to work, from my experimentation. But the values are then "frozen" in the machine. Even going to a different script does not enable those variable names to be changed. I'm not too impressed with Octave right now!