Declaring an anonymous function for solving an odes system

13 views (last 30 days)
Hi there,
I am trying to integrate a system of ODES, I am being partially successful defining the function in another file, but when trying to do that using an anonymus function, it is not working, I am not sure what is wrong.
So, the working version is:
in F.m file:
function x = F(t,theta)
N = 10; lambda = 1; k = 4;
omega = lambda*tan(pi*(rand(1,N)-0.5));
x = omega' + k/N*sum(sin(repmat(theta,[1 N])' - repmat(theta,[1 N])),2);
and then:
N = 10;
[t,x] = ode45(@F,[0,10],ones(N,1));
And the one Im trying to make work is:
N = 10; lambda = 1; k = 4;
omega = lambda*tan(pi*(rand(1,N)-0.5));
FK = @(t,theta) omega' + k/N*sum(sin(repmat(theta,[1 N])' - repmat(theta,[1 N])),2) ;
[t,x] = ode45(FK,[0,10],ones(N,1));
Thanks in advance
  2 Comments
Adam
Adam on 29 Sep 2015
What do you mean by 'not working'? Is there an error message or are you just getting incorrect results?
gringer45
gringer45 on 29 Sep 2015
Incorrect results, totally different to the results I get with the other method. No matlab errors though.

Sign in to comment.

Accepted Answer

gringer45
gringer45 on 29 Sep 2015
Ok, I just figuered it out.
The omega function was being defined only once in the second case, so it was being used the same for each run instead of a new one.
FK = @(t,theta) lambda*tan(pi*(rand(N,1)-0.5)) + k/N*sum(sin(repmat(theta,[1 N])' - repmat(theta,[1 N])),2) ;
instead of
omega = lambda*tan(pi*(rand(1,N)-0.5));
FK = @(t,theta) omega' + k/N*sum(sin(repmat(theta,[1 N])' - repmat(theta,[1 N])),2) ;
cheers
  1 Comment
Walter Roberson
Walter Roberson on 29 Sep 2015
No, the other way around. When you use rand() or randn() in your ODE function, you make the function discontinuous and non-repeatable. This would drive the integrator to smaller and smaller step sizes unless you are using a fixed-step solver. If you must have randomness, you should define the random contributions ahead of time, possibly with some method of interpolating at interior times.

Sign in to comment.

More Answers (0)

Products

Community Treasure Hunt

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

Start Hunting!