0

I'm new in Matlab, so sorry for asking a banal and probably already answered question, but I wasn't able to solve it.

I have this ODE to solve

g''(x)+f(x)g'(x)

with g(0)=0 and g'(0)=0.5

f(x) is a known vector, so I tried

xspan = linspace(0,10,100);
y0 = [0,0.5];
%f=f(xspan), known function, for example f=1./xspan;
[x,y] = ode45('odefun',xspan,y0);

and for odefun 

function dy=odefun(x,y)

global f  

dy(1) = y(2);
dy(2) = -f.*y(2);

dy    = dy(:);
end

Of course it doesn't work. How can I pass the values of vector f to the solver?

EDIT

Thanks for help, but this is not working for me. Maybe I have an old Matlab version (7.11.0 (R2010B))? I tried to do as you told, but this is what I get in command window

???  In an assignment  A(I) = B, the number of elements in B and
 I must be the same.

Error in ==> odefun at 3
dy(2) = -f.*y(2);

Error in ==> @(x,y)odefun(x,y,f)


Error in ==> odearguments at 109
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ==> ode45 at 173
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Error in ==> trying at 7
[x,y] = ode45(fun,xspan,y0);

This is the code trying.m

xspan=linspace(0,10,100);
y0=[0 0.5];
f=1./xspan;
[x,y] = ode45(@(x,y)odefun(x,y,f),xspan,y0);

and the odefun

function dy=odefun(x,y,f)
dy(1) = y(2);
dy(2) = -f.*y(2);
dy    = dy(:);
end
user149839
  • 3
  • 2
  • Re: your edit (a new question): You still have to satisfy all of the usual rules for writing ODE functions. It looks like you're passing in a vector for `f` and then trying to set it equal to a single element of `dy`: `dy(2) = -f.*y(2);`. `dy(2) = -f(1).*y(2);` or `dy(2) = -f(end).*y(2);` would work. Are you possibly trying to extract the step size, `dt` that `ode45` is using on each time step to use it in `odefun`? If so, this approach won't work. – horchler May 16 '14 at 18:05
  • I didn't think about that, f is a vector and I can't do that of course. Sorry for the mistake. By the way I don't want to do f(1)*y(2) or f(end), I need f to be evaluated at every step used in the ode45 solver. So, should I find the step size, as suggested, and then interpolate the f at those steps? How could I do that? Really thank you for your help! – user149839 May 16 '14 at 20:52
  • There is no way to get the step sizes that `ode45` is using internally while it runs (`xspan` is only the set of output points not necessarily the true steps taken by the solver – [my answer here](http://stackoverflow.com/a/21861868/2278029) might help a bit). You can get the independent variable `x`, but not the step size. An ODE shouldn't be a function of step size to begin with, only the current independent and dependent variables and any parameters. You may be trying something numerically and mathematically inadvisable. – horchler May 16 '14 at 21:03
  • In the truth I don't need to know the step sizes of the ode45, I need to solve the differential equation g''(x)+f(x)g'(x)=0, and I didnt know how the pass the values of the f(x) (that I evaluated before). I understand that the xspan vector does not represent the true steps used by the ode45, but I needed the f(x) to be evaluated at those steps to solve the equation above, for this I thought about interpolation. So, is there no way I can do this? If I can get the independent variable x, so could I use interpolation to get f(x) values for each 'x' of the ode? If not, is there no way to solve it? – user149839 May 17 '14 at 13:17
  • The first input to your `odefun` is `x`. – horchler May 17 '14 at 15:43
  • If you could help me with the code I would apprecciate, because as I told you I'm not very used with programming. So how can I get the values of f(x) knowing f(xspan)? – user149839 May 17 '14 at 18:53
  • This is really a separate question from the original one that you asked. Also, you haven't defined your function `f` as a function. If f(x) is a constant function that doesn't actually depend `x` then you can pass in a scalar or vector. Otherwise you need to pass in a function (i.e., a function handle) that takes `x` as an argument. You need to define exactly what `f` is and ask a new specific question related to this. Provide any error messages that you get for what you're trying, etc. – horchler May 17 '14 at 20:39
  • f(x) is a function get from the solution fo another ode, so I know it at each xspan, even if I don't know its analitical expression. If you want I asked a related question here [link](http://stackoverflow.com/questions/23737271/how-to-get-values-by-interpolation-of-a-function-fx-not-knowing-its-analit) thanks for help – user149839 May 19 '14 at 12:07
  • It was very easy. In the way I put in the answer the code works, but why if `f` is a vector, and I set it equal to a single element, as we observed above? It works and it's good for me , but I can't get it :P – user149839 May 21 '14 at 21:46

2 Answers2

1

It looks like you're using global variables incorrectly (f needs to be declared global in your main function as well), but you shouldn't be using globals to begin with. You should rarely, if ever, use globals. You should also specify your integration function, odefun, as an anonymous function rather than a string in modern Matlab. Here is how you can modify your code to fix these issues:

...
f = ...
[x,y] = ode45(@(x,y)odefun(x,y,f),xspan,y0);

And your integration function, which can be a sub-function or in a separate file, just needs to take an extra argument for the parameter f:

function dy=odefun(x,y,f)
dy(1) = y(2);
dy(2) = -f.*y(2);
dy    = dy(:);

The methods should also be faster than using strings and globals. The code @(x,y)odefun(x,y,f) creates an anonymous function of x and y (ode45 expects a function that takes exactly two inputs) that captures the value of f from the current scope. This is also referred to as a closure in computer science. You can equivalently create a function handle and pass that in:

...
f = ...
fun = @(x,y)odefun(x,y,f)
[x,y] = ode45(fun,xspan,y0);

See also this article on parametrizing functions from The MathWorks for more information on this general technique.

horchler
  • 18,384
  • 4
  • 37
  • 73
  • Thanks for interest, but I've still got problems. I edited, the question, if you could I would appreciate your help – user149839 May 16 '14 at 13:21
0

I solved it, it was not so difficult, I just needed to interpolate the f at each value of the independent variable x. In this way the code works

function dy=odefun(x,y,f,xspan)
f=interp1(xspan,f,x);
dy(1) = y(2);
dy(2) = -f.*y(2);
dy    = dy(:);
end
user149839
  • 3
  • 2