Plot dB/dt vs. time by ODE45 (with dB/dt on the y axis and time on the x axis)

3 ビュー (過去 30 日間)
S H
S H 2019 年 12 月 14 日
コメント済み: S H 2019 年 12 月 17 日
Hi,
I have plotted the solution for two odes against time.
The two odes are given in the following.
dA/dt=-k1*(exp(-E1/RT))*A
dB/dt=k1*(exp(-E1/RT))*A-k2*(exp(-E2/RT))*B
T=25+5t
E1=60
k1=1e-2
E2=100
k2=1e-2
R=8.314
What I would like to do is also plot dB/dt vs. time (differential equations themselves).
The code I have now is not wokring properly. I was wondering if I make any mistake with the codes?
function diffeqs=odes(t,var)
format long
T=var(1);
A=var(2);
B=var(3);
E1=60;
k1=1e-2;
E2=100;
k2=1e-2;
R=8.314;
if t < 20
diffeqs(1,1)=5;%dT/dt
else
diffeqs(1,1)=0;%dT/dt
end
diffeqs(2,1)=-k1*exp(-E1/(R*T))*var(2);%dA/dt
diffeqs(3,1)=k1*exp(-E1/(R*T))*var(2)-k2*exp(-E2/(R*T))*var(3);%dB/dt
end
Plotting
clc;close all; clear all;
timerange=linspace(0,100,50);
ICs= [25,10,0];
[tsol,varsol]=ode45(@(t,var) odes(t,var), timerange, ICs);
pvar = gradient(varsol(:,3));
dt = gradient(tsol);
figure
yyaxis right
plot(tsol,varsol(:,1),'r','linewidth',2);
xlabel('Time, min');
ylabel('Temperature');
hold on
yyaxis left
plot(tsol,varsol(:,3),'b','linewidth',2);
xlabel('time');
ylabel('B(t)');
grid
legend('B' );
figure
hq = quiver(tsol, varsol(:,3), dt, pvar, 6);
set(hq,'Color','k','linewidth',2)
xlabel('time');
ylabel('dB/dt');
grid
col_header1={'Time','B'}
xlswrite('Book2.xlsx',[tsol(:),varsol(:,3)],'sheet1','A2');
xlswrite('Book2.xlsx',col_header1,'sheet1','A1');
Figure 1 is the plot I got for B(t) vs. time.
Figure 2 is dB/dt vs. time (which is incorrect result)
If I extract the datapoints from Figure 1 and replot it in excel. I got Figure 3 which is the same as Figure 1.
After that, I was using these datapoints to plot dB/dt vs. time then got Figure 4 which is very different from Figure 2.
Therefore, I was wondering if there is anything wrong with my codes? It would be greatly appreciated if anyone can help. Thank you in advance!

採用された回答

darova
darova 2019 年 12 月 14 日
編集済み: darova 2019 年 12 月 14 日
  • What I would like to do is also plot dB/dt vs. time
pvar = gradient(varsol(:,3));
dt = gradient(tsol);
dB = pvar./dt;
plot(tsol,dB)

その他の回答 (0 件)

カテゴリ

Help Center および File ExchangeAnnotations についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by