ode within fsolve routine

Hello Matlab-community,

I have a problem with a physical model of a membrane. I have to integrate a set of derivatives from the back end (omega) to the front end (alpha) of the membrane, compare the calculated values at the front end with the actual values and modify the initial condition (at the back end) with a solver until the calculated vales are the same as the measured values.

My program has the following structure:

 % main program

% final values
yyalpha=[1;2;3;4];

%estimated initial conditions
yyomega0=0.5*yyalpha;

%solver for initial conditions
yyomega=fsolve(targetfun,yyomega0)
disp(x,yyomega);

The main program calls the target function. It is:

%targetfun
function tf=targetfun(yyomega,yyomega0,yyalpha)
    [x,yy]=ode45(yydot,0:0.1:1,yyomega);
    tf=yy(:,11)-yyalpha;
end

and the target function calls the derivative:

%derivative
function dyy=yydot
    dyy=yyalpha+yy;
end

The problem is that the variables are unknown in the functions.

Error in yydot (line 3)
    dyy=yyalpha+yy;

Error in targetfun (line 4)
    [x,yy]=ode45(yydot,0:0.1:1,yyomega);

Error in simpleexample (line 10)
yyomega=fsolve(targetfun,yyomega0)

I changed the code several times so that all functions handle through the variables but that didn't help. there are always undefined variables. When I handle through all the variables the solver variable (yyomega) is unknown. I also used function handles (@) but that didn't help either.

How can I change the code so the calculation is carried out? Is it possible to nest an ODE into a fsolve-routine?

Answer 1

[Walter Roberson]

% main program

% final values
yyalpha = [1;2;3;4];

%estimated initial conditions
yyomega0 = 0.5*yyalpha;

%solver for initial conditions
yyomega = fsolve(@(yyomega) targetfun(yyomega, yyalpha), yyomega0);
disp(yyomega);       %x is not easy to return from targetfun

%targetfun
function tf = targetfun(yyomega, yyalpha)
    [x, yy] = ode45(@(x, yy) yydot(x, yy, yyalpha), 0:0.1:1, yyomega);
    tf = yy(end, :).' - yyalpha;
end

%derivative
function dyy = yydot(~, yy, yyalpha)
  dyy = yyalpha + yy;
end

Answer 2

[JK]

ok, thanks for this. It seems I have to dig a little deeper into the arcane art of using function handles.

My deductions are:

1. all functions must be called with the full set of variables which should be passed from one function layer to another. Variables which are defined locally at function layer 1 are not valid in function layer 2 (subfunction)
2. the variables which should be given back by the function must be marked with the function handle (@) and stated before the function name when it is called.
3. the ~ is just for show off. it also works if the ~ is replaced with the variable name (here x)
4. ode45 seems to give back a row-vector and not a column-vector

once again thanks for the help.

[Walter Roberson]

1) Subfunctions can access variables defined in their parent function that were initialized in the code before the subfunction was defined. However, your code does not nest functions, so you have no subfunctions.

   function outer
     some_variable = 17.982;  %must be defined before inner function is declared

     function y = inner
       y = some_variable.^2;   %inner function can access variables in outer function that were assigned to before inner was defined
     end

     some_other_variable = 192.318;
     yet_other_variable = inner();   %this is a continuation of function outer

   end

function inner would not have access to some_other_variable even though some_other_variable is assigned to before the call to the inner function: only variables that have been assigned to before the inner function body appears can be shared with the inner function.

But as I said, you did not nest functions, so no sharing of variables is applicable to what you had posted.

2) "The variables which should be given back by the function must be marked with the function handle (@) and stated before the function name when it is called."

No, that is not the case. Anonymous function syntax is

@(variable1, variable2 ...) expression_to_execute(....)

where variable1, variable2 ... are the dummy names of the positional parameters that the anonymous function will be called with. The expression_to_execute that is invoked can use any of those dummy parameter names in any order, and can also refer to constants and other variables. If other variables are referred to, then those other variables have to have been assigned to already (or be the names of parameters that were passed in to the current function), and a copy of the value of those variables is taken as of the time the anonymous function is formed, and those copies are bound together with the anonymous function. For example,

yyalpha = 1:10;
fun = @(x) norm(x-yyalpha);
yyalpha = 1492;
fzero(fun, [0 48])

Here x is the dummy parameter name referring to the first positional parameter that fun will be invoked with. The additional variable yyalpha is referred to, so its value at the time, 1:10, will be copied and bound together with the code branch, creating a workspace with variable yyalpha and value 1:10. The change to yyalpha on the example code line after that has no effect on fun: the value of yyalpha was already copied. Anonymous functions do not reference variables by name such that later changes to variables with that name might have impact: anonymous functions copy the value at the time of definition. So when the anonymous function fun is invoked in the fzero context, the yyalpha value for the purposes of execution will be the 1:10 value.

This has nothing to do with "variables given back". It would have been completely valid to have used

solution_yyomega = fsolve(@(yyomega) targetfun(yyomega, yyalpha), yyomega0);
disp(solution_yyomega);       %x is not easy to return from targetfun

The dummy variable names used in defining the anonymous functions do not need to be used anywhere else.

3. With regards to the tilde being for show-off: the optimizer can do a better job when the tilde is used in place of a variable whose value will not be used. I originally wrote the code using the output name x there, but the MATLAB editor reminded me that it was better to use ~ if I was not going to use the result later.

4. "ode45 seems to give back a row-vector and not a column-vector"

Each of the output variables from the column output by the invoked function is made into a column. When there are multiple variables, the result will almost always be a 2D array, not a vector (the event function would have to have fired on the first time to get less than two rows of outputs.) There will be one row for each output time. You happened to be interested only in the last of those, yy(end,:) .

I agree there is a slight inconsistency between requiring that the invoked function always return a column vector, but then turning around and making that into a row.

Note along the way, though: the outputs from ode45 are not typically just the results of the individual invocations or even just the integral of the results of the individual invocations. Instead, ode45 makes several calls to determine conditions "near" the point of interest, and uses Range-Kutta(4,5) to predict the output that it uses.

Answer 3

[JK]

wow, if you have only spend half the time for writing as I have been spending for reading and trying to understand you have made quite an effort. I really appreciate your time.

the function in your example should read

fun = @(x) norm(x-yyalpha)-50

otherwise it won't get negative and fzero cannot find the root. But I get your point. I have programmed a little in my youth - C64 BASIC, turbo pascal and assembler, but Matlab is so much more complicated. I am currently using Matlab for Dummies as a learning book but I am not very satisfied. It focuses for ever on Hello World-stuff but does not explain the essentials like why it is important to use these anonymous functions. Can you recommend a good book which focuses on the middle part between the very simple and very complicated matter?

[Walter Roberson]

You are right, that particular example with norm should subtract something off to have a zero crossing. I was not paying attention to whether the function had a root; I was just tossing code together to illustrate the concepts.

I do not have any particular books to recommend. A while ago people put togther http://www.mathworks.com/matlabcentral/answers/8026-best-way-s-to-master-matlab .