1

I'm solving an ODE like

y''(t) + a y(t) + b = 0

with a Matlab's ode45 solver. It iterates until a solution was found, which satisfies the error criteria. I want to read the previous value of y in each step for a comparison.

What would be the best way to do this?

My main file is

[t,y] = ode45(@odefnc,tspan,[0 0]',options);

and the odefnc is

function dx=odefnc(t,x)

...

if history(end)<13
    dx=[x(2),-a*x(1)-b]'
else
    dx=[x(2),-c*x(1)-d]'
end


if flag==1
   history(end+1)=dx;
end
horchler
  • 18,384
  • 4
  • 37
  • 73
Caniko
  • 867
  • 2
  • 11
  • 28

1 Answers1

1

You're trying to change the ODEs from within the integration function by adding large discontinuities (the if statements). This is bad practice and can lead to all sorts of issues. Don't do it. Instead you need to call ode45 multiple times and change the the ODE function – or the parameters you give it. See my answer here for further details and some examples.

You also seem to be trying to record the "history" of the evaluated solution. You can't do this. You cannot rely on the order that your ODE function gets called by ode45 or how many times it gets called per time step (maybe many times in the case a failed step). What you seem to be trying to do just is not possible with ode45 and the other ODE Suite functions. However, it's likely that whatever want to do can be accomplished in some better, more natural way. Or you may need to implement your own integration scheme – but this is rarely necessary.

If all you're trying to do by recording the value of the previous step, is determine when to switch between the parameter a and the parameter c, then you need a way of determining that point accurately (your if statements, even if they work, will not be accurate). If the switch occurs at a particular point in time, then it's easy: just perform two integrations over the necessary time spans as shown at the link above. If the switch occurs at a condition that depends on the state x, then you'll need to learn about event functions. You can read about those in the help and documentation, this article from the MathWorks, and also my answers to this question and this question.

Community
  • 1
  • 1
horchler
  • 18,384
  • 4
  • 37
  • 73
  • Horchler, thank you for your editing and detailed answer on my post! The fact is, I am modelling a mechanical system with both slipping and adherence, that means also their transitions in each other so that the large discontiniuties aren't avoidable for me. The system shall work within any other global multibody system embedded which causes that I won't be able to predict the transitions (system changes). I already gained consistent results embedding my tool into a commercial mbs-solver as global system solver. Now I want to write my own MBS Solver in Matlab for a better performance... – Caniko Aug 07 '13 at 09:51
  • ...This will solve a very simple global system (saying a box with mass) However even in this case I will append my discontinous mechanical system as an external routine solving itself undependent but returning discont. results as expecting. I mean, theoretical this should be possible.. – Caniko Aug 07 '13 at 09:56
  • 1
    @Caniko: Everything I wrote still stands. Consistency is not accuracy in numerical integration. You will need to write event detection functions to find the points at which the system slips, stick, and is totally free. It's tricky. I really suggest you read some on this problem rather than implementing it haphazardly. Or you might look into Matlab's [SimScape](http://www.mathworks.com/products/simscape/) environment. In any case `ode45` is may well be the wrong solver for this. You should look into [`ode15`](http://www.mathworks.com/help/matlab/ref/ode15s.html). – horchler Aug 07 '13 at 14:24
  • thank you for your useful answers! After researching according to your tipps I think I know how to approach my problem! It's more or less about starting a new problem(ode) with new initial conditions by detection of an event change. You're absolutely right. I also faced the mentioned matters implementing my problem. – Caniko Aug 08 '13 at 09:16