- What I would like to do is also plot dB/dt vs. time
Plot dB/dt vs. time by ODE45 (with dB/dt on the y axis and time on the x axis)
3 ビュー (過去 30 日間)
古いコメントを表示
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 件のコメント
採用された回答
その他の回答 (0 件)
参考
カテゴリ
Help Center および File Exchange で Annotations についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!