How can i plot dy/dx function?

24 次查看（过去 30 天）
nur farizan munajat
编辑： Torsten 2022-10-8
Ran in:
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

请先登录，再回答此问题。

采纳的回答

Davide Masiello
Davide Masiello 2022-10-6
编辑：Davide Masiello 2022-10-6
Ran in:
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.

更多回答（1 个）

Torsten
Torsten 2022-10-6
Ran in:
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 个评论
nur farizan munajat
nur farizan munajat 2022-10-8
Thank you for your answer. Now i expand my function dydx. This funtion is actually based on the experimental data. it suppose to get more than one peak. but i still got only one peak
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;
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);
dydT =(Q1+Q2+Q3+Q4)/5;
end
Torsten
Torsten 2022-10-8
编辑：Torsten 2022-10-8
Ran in:
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

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 MATLAB 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by