Declaring an anonymous function for solving an odes system

13 visualizzazioni (ultimi 30 giorni)
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 Commenti
Adam
Adam il 29 Set 2015
What do you mean by 'not working'? Is there an error message or are you just getting incorrect results?
gringer45
gringer45 il 29 Set 2015
Incorrect results, totally different to the results I get with the other method. No matlab errors though.

Accedi per commentare.

Risposta accettata

gringer45
gringer45 il 29 Set 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 Commento
Walter Roberson
Walter Roberson il 29 Set 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.

Accedi per commentare.

Più risposte (0)

Prodotti

Community Treasure Hunt

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

Start Hunting!

Translated by