0

I want to solve an ode with a time dependent parameter. cA should be 10000 if t is >=10 and <=11 else it should have the value of 0. cA is then used in an differential equation to calculate cB. See the code:

function dcB = myode(t,y)
cB=y(1,:);

if t>=10 && t<=11
    cA=10000
else
    cA=0
end
dcB=(cA-cB)*100/1750;

[t,y]=ode45(@myode,[tdown tup],0);

Fallowing problems show up:

  • if I print cA it has not the correct values at the specified times.
  • if tup is e.g. 20 cB has values, if tup is e.g. 100 cB is zero.
Josh Crozier
  • 233,099
  • 56
  • 391
  • 304
sejo
  • 1
  • 4
  • If I modify your `myode` so that it prints both `t` and `cA` each time the function is called, then call it using `[t,y]=ode45(@myode,[0 20],0);plot(t,y);` then the numbers printed to the MATLAB Command Window appear to be correct, and the plot shows the expected discontinuity at 10 and 11 seconds. – Phil Goddard Jul 23 '14 at 15:15
  • What you're trying to do inside of your `myode` integration function with the `if` statement is a bad idea. User-supplied ODE functions should not have discontinues. It will only lead to inefficient computation and inaccuracies. See my answers to [this question](http://stackoverflow.com/a/17370733/2278029) and [this one](http://stackoverflow.com/q/21500204/2278029) for the proper way to accomplish what you're trying to do. If you want a "time-dependent parameter", it should vary continuously. – horchler Jul 23 '14 at 18:34
  • @PhilGoddard yes with tup=20 it works but if you try with 100 the result is zero... – sejo Jul 25 '14 at 08:50
  • The maximum step size is most likely too large, and the solver is stepping over your critical points. Try `options = odeset('MaxStep',0.01);` or some other suitable step size and then `[t,y]=ode45(@myode,[0 100],0,options);plot(t,y);` – Phil Goddard Jul 25 '14 at 14:03

1 Answers1

0

As horchler suggested i solved the problem according to this answer: link

function dcB = myode(t,y,cA)
cB=y(1,:);
dcB=(cA-cB)*100/1750;

and solve with

yO=0;

% t=0 till t=10
cA=0;
tspan=[0 10];
[T,Y]=ode45(@(t,y)myode(t,y,cA),tspan,yO,cA);

% t=10 till t=11
cA=10000;
tspan=[0 1];
[t,y]=ode45(@(t,y)myode(t,y,cA),tspan+T(end),Y(end,:));
T = [T;t(2:end)];
Y = [Y;y(2:end,:)];

% t=11 till t=100
cA=0;
tspan=[0 89];
[t,y]=ode45(@(t,y)myode(t,y,cA),tspan+T(end),Y(end,:));
T = [T;t(2:end)];
Y = [Y;y(2:end,:)];

plot(T,Y);
Community
  • 1
  • 1
sejo
  • 1
  • 4
  • Glad to see you solved your own question and thank you for posting the solution. If any of the answers at the link you provided were helpful, please [vote them up](http://stackoverflow.com/help/why-vote). – horchler Jul 28 '14 at 14:50