MATLAB Answers

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

3 views (last 30 days)
S H
S H on 14 Dec 2019
Commented: S H on 17 Dec 2019
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!

  0 Comments

Sign in to comment.

Accepted Answer

darova
darova on 14 Dec 2019
Edited: darova on 14 Dec 2019
  • 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)

More Answers (0)

Sign in to answer this question.


Translated by