Declaring an anonymous function for solving an odes system

1 次查看（过去 30 天）
gringer45
gringer45 2015-9-29
评论： Walter Roberson 2015-9-29
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 个评论
Adam
Adam 2015-9-29
What do you mean by 'not working'? Is there an error message or are you just getting incorrect results?
gringer45
gringer45 2015-9-29
Incorrect results, totally different to the results I get with the other method. No matlab errors though.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

gringer45
gringer45 2015-9-29
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 个评论
Walter Roberson
Walter Roberson 2015-9-29
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.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by