How can i plot dy/dx function?
41 views (last 30 days)
Show older comments
nur farizan munajat
on 6 Oct 2022
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
0 Comments
Accepted Answer
Davide Masiello
on 6 Oct 2022
Edited: Davide Masiello
on 6 Oct 2022
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.
0 Comments
More Answers (1)
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
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)
hold off
figure(2)
hold on
plot(x1,dy1)
plot(x2,dy2)
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
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!