How can i plot dy/dx function?
42 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
nur farizan munajat
am 6 Okt. 2022
Bearbeitet: Torsten
am 8 Okt. 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 Kommentare
Akzeptierte Antwort
Davide Masiello
am 6 Okt. 2022
Bearbeitet: Davide Masiello
am 6 Okt. 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 Kommentare
Weitere Antworten (1)
Torsten
am 6 Okt. 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 Kommentare
Torsten
am 8 Okt. 2022
Bearbeitet: Torsten
am 8 Okt. 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
Siehe auch
Kategorien
Mehr zu Ordinary Differential Equations finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!