0
votes

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.
1
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
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 and this one for the proper way to accomplish what you're trying to do. If you want a "time-dependent parameter", it should vary continuously.horchler
@PhilGoddard yes with tup=20 it works but if you try with 100 the result is zero...sejo
@horchler Thanks for your answer. I will try it that way!sejo
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

1 Answers

0
votes

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);