How can i plot dy/dx function?

29 views (last 30 days)
i have below code. i can get the solution for the derivative equation and able to plot y vs x. But how can i plot dydx vs x?
clear
xspan = (30:0.1:900);
y0 = 0;
[x,y] = ode23(@(x,y) odefun(x,y), xspan, y0);
plot(x,y)
function dydx = odefun(x,y)
R=8.314;
Q1=0.29*7.54E19*exp(-209201/(R*(x+273)))*(1-y)^3;
Q2=0.34*5.07E16*exp(-195508/(R*(x+273)))*(1-y)^2;
dydx =(Q1+Q2)/5;
end

Accepted Answer

Davide Masiello
Davide Masiello on 6 Oct 2022
Edited: Davide Masiello on 6 Oct 2022
You can use deval in combination with solution structure.
xspan = (30:0.1:900);
y0 = 0;
opts = odeset('AbsTol',1e-8,'RelTol',1e-6);
sol = ode15s(@(x,y) odefun(x,y), xspan, y0, opts);
[~,yp] = deval(sol,sol.x);
plot(sol.x,sol.y)
ylabel('y')
yyaxis right
plot(sol.x,yp)
ylabel('dy/dx')
function dydx = odefun(x,y)
R=8.314;
Q1=0.29*7.54E19*exp(-209201/(R*(x+273)))*(1-y)^3;
Q2=0.34*5.07E16*exp(-195508/(R*(x+273)))*(1-y)^2;
dydx =(Q1+Q2)/5;
end
PS. I changed solver and adjusted tolerances for a better result, you might want to consider doing the same.

More Answers (1)

Torsten
Torsten on 6 Oct 2022
xspan = (30:0.1:900);
y0 = 0;
options = odeset('RelTol',1e-12,'AbsTol',1e-12);
[x,y] = ode15s(@(x,y) odefun(x,y), xspan, y0, options);
dy = zeros(size(y));
for i = 1:numel(x)
dy(i) = odefun(x(i),y(i));
end
%hold on
%plot(x,y)
plot(x,dy)
%hold off
function dydx = odefun(x,y)
R=8.314;
Q1=0.29*7.54E19*exp(-209201/(R*(x+273)))*(1-y)^3;
Q2=0.34*5.07E16*exp(-195508/(R*(x+273)))*(1-y)^2;
dydx =(Q1+Q2)/5;
end
  2 Comments
Torsten
Torsten on 8 Oct 2022
Edited: Torsten on 8 Oct 2022
The solution curve for y shows that there should only be one peak ...
Obviously, one reaction dominates the other three massively:
xspan = (30:0.1:900);
y0 = 0;
options = odeset('RelTol',1e-12,'AbsTol',1e-12);
info = 1;
[x1,y1] = ode15s(@(x,y) odefun(x,y,info), xspan, y0, options);
info = 2;
[x2,y2] = ode15s(@(x,y) odefun(x,y,info), xspan, y0, options);
dy1 = zeros(size(y1));
for i = 1:numel(x1)
dy1(i) = odefun(x1(i),y1(i),1);
end
dy2 = zeros(size(y2));
for i = 1:numel(x2)
dy2(i) = odefun(x2(i),y2(i),1);
end
figure(1)
hold on
plot(x1,y1)
plot(x2,y2)
Warning: Imaginary parts of complex X and/or Y arguments ignored.
hold off
figure(2)
hold on
plot(x1,dy1)
plot(x2,dy2)
Warning: Imaginary parts of complex X and/or Y arguments ignored.
hold off
function dydT = odefun(x,y, info)
R=8.314;
Q1=0.29*7.54E19*exp(-209201/(R*(x+273)))*(1-y)^3;
Q2=0.34*5.07E16*exp(-195508/(R*(x+273)))*(1-y)^2;
Q3=0.12*1.48E07*exp(-136978/(R*(x+273)))*2*(1-y)*sqrt(-log(1-y));
Q4=0.09*6.88E07*exp(-142875/(R*(x+273)))*3*(1-y)*(-log(1-y))^(2/3);
if info == 1
dydT = Q1/5;
else
dydT = (Q1+Q2+Q3+Q4)/5;
end
end

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!

Translated by