9 views (last 30 days)

Hi!

I am using a function with a subfunction to produce a graph of a piecewise smooth variable in an ODE. Smoothness is lost when t=t_1. Code for this is shown at the bottom of the question.

I'd like code to produce a 3d graph showing how different solutions look as t_1 varies. This would have the three axes t, t_1 and xa(:,1). Currently what I am stuck on is I am not sure how I could include t_1 as a a variable that is fed into ydot without creating confusing when ode45 calls it.

function L2_ode45

%a function to call ode45 and draw the graph

[t,xa] = ode45(@fu,[0 500],[0 0]);

plot(t,xa(:,1))

grid on

end

function ydot = fu(t,x)

%the function passed to ode45

%give values to parameters

t_1 = 200

k_a1 = 10^(-4);

k_a2 = 2*10^(-2);

C = 1;

R = 10;

k_d1 = 10^(-3);

k_d2 = 10^(-3);

%set up the switch

if t < t_1

z = C*(R-x(1)-x(2));

else

z = 0;

end

%set up the ROC function that will be output to ode45

ydot = [ k_a1*z-k_d1*x(1);...

k_a2*z-k_d2*x(2);];

end

Star Strider
on 15 Jan 2020

Edited: Star Strider
on 15 Jan 2020

I am not certain what you want to do.

Try this:

function L2_ode45

tspan = linspace(0, 500, 50);

t1 = 100:100:500;

for k = 1:numel(t1)

[t{k},xa{k}] = ode45(@(t,x)fu(t,x,t1(k)),tspan,[0 0]);

end

q1 = size(t)

tp = cell2mat(t);

xap = cell2mat(xa);

xapm = xap(:,1:2:end);

figure

mesh(t1, tspan, xapm)

grid on

xlabel('t_1')

ylabel('t')

end

function ydot = fu(t,x,t_1)

%give values to parameters

% t_1 = 200

k_a1 = 10^(-4);

k_a2 = 2*10^(-2);

C = 1;

R = 10;

k_d1 = 10^(-3);

k_d2 = 10^(-3);

%set up the switch

if t < t_1

z = C*(R-x(1)-x(2));

else

z = 0;

end

%set up the ROC function that will be output to ode45

ydot = [ k_a1*z-k_d1*x(1);...

k_a2*z-k_d2*x(2);];

end

If you want to interrupt the integration to substitute some new value, see the documentation on ODE Event Location.

EDIT —

I forgot to post the plot!

Here it is —

Sign in to comment.

Sign in to answer this question.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.