This problem involves 3 time dependent ODEs that are to be solved graphically with the input of a time dependent function.

1 次查看（过去 30 天）
James
James 2011-4-1
Dear All,
I briefly came across the following problem in Matlab:
Solve the following ODEs:
dA/dt = k2B - k1A; dB/dt = -k2B + k1A - k3B; dC/dt = k3B;
where: k1=A1*exp(-Ea1/(R*(T+273))); k2=A2*exp(-Ea2/(R*(T+273))); k3=A3*exp(-Ea3/(R*(T+273))); and: R=8.3145 Ea1=50000; Ea2=Ea1; Ea3=25000; A1=8.2e10; A2=8.2e9; A3=2409;
and we also have the time dependent temperature:
T(t)= 20 + 8t; for 0 < t < 10.
What is the best way to solve this problem? I believe a separate function along with an M file is required. Also, I usually use ODE45 for solving such ODEs. To input the different T into k1, k2, k3, is a for loop required or is an array better?
Any help would be greatly appreciated. Thanks in advance,
James Champion

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

回答（3 个）

Matt Tearle
Matt Tearle 2011-4-1
When using ode45, it will evaluate your function at a given set of values for the independent and dependent variables, so you don't need a loop. That is, the looping is done in ode45, not your function.
So write a function file
function dy = myodefun(t,y)
Temp = 20 + 8*t;
k1 = ...
[etc]
dy = zeros(3,1);
dy(1) = k2*y(2) - k1*y(1);
[etc]
then pass that to ode45. Note the use of the vectors y and dy to represent the three dependent variables A, B, and C, and their derivatives.
James
James 2011-4-2
Hi Matt,
Thanks for the response. If I wrote my function like this:
function xdot = odefun(t,x) T = 20 + 8*t;
R=8.3145;
Ea1=50000; Ea2=Ea1; Ea3=25000;
A1=8.2e10; A2=8.2e9; A3=2409;
xdot=zeros(3,1);
xdot(1)=(-k1*x(1))+(k2*x(2));
xdot(2)=(k1*x(1))-(k2*x(2))-(k3*x(2));
xdot(3)=(k3*x(2));
end
and an ODE45/mfile like this:
global k1 k2 k3
t0=0;tf=10;
x0=[10 100 0];
[t,x]=ode45('odefun',t0,tf,x0);
plot(t,x)
title ('odefun')
would this solve the original time dependent problem? It is difficult for me to test as I currently don't have access to Matlab.
Thanks,
James
  1 个评论
Matt Tearle
Matt Tearle 2011-4-3
Mostly looks ok, but I don't understand why you're using global for k1, k2, k3 (or where those are being calculated). Why not just calculate them inside your function? You've calculated T, and set all the constants (R, Ea1, A1, etc), so you have everything you need to get the ks as well.

请先登录，再进行评论。

Walter Roberson
Walter Roberson 2011-4-3
James sent the problem to me in email as well. Working with the problem with another problem, it does not appear that there are useful solutions due to the lack of initial conditions (unless his A1 means A(1))
  3 个评论
显示 1更早的评论
Walter Roberson
Walter Roberson 2011-4-4
With the initial conditions, the tool I am using is happy to give a try at numeric solutions. The fastest so far is Livermore Stiff with Jacobian, or perhaps Modified Extended BDF Explicit.
The fixed-step methods produce some pretty crazy answers.
Gear single step extrapolation is really the only sort-of-reasonable solution that has a notable distinguishing feature, having a deep valley near 5.8 that there is not even a hint of on anything else.
Walter Roberson
Walter Roberson 2011-4-4
The symbolic solution has been calculating for a rather large number of hours. I would not count on a reasonable symbolic solution being forthcoming.
You might wish to consider the following equivalence, which is easier for some of the tools to work with:
diff(A(t),t) = k2(t)*B(t) - k1(t)*A(t);
diff(B(t),t) = -diff(A(t),t) - diff(C,t);
diff(C(t),t) = k3(t)*B(t);
k1(t) = A1*exp(-Ea1/(R*(T(t)+273)));
k2(t) = A2*exp(-Ea2/(R*(T(t)+273)));
k3(t) = A3*exp(-Ea3/(R*(T(t)+273)));
T(t) = 20 + t*8;
Along with the parameters.
As you are working with floating point numbers, the additional hint that the initial values for diff(B(t),t) are not just "close to" those of diff(A(t),t) but are _identical_ to them, can help reduce the complexity of the system.
Don't even bother considering any of the fixed step-size numeric routines, though -- they all go wild somewhere near 6.
Taylor expansion of order 6 is fairly practical but gives obviously poor answers; expansion of order 11 takes notably more computation by still gives obviously poor answers.
Basically if you use any of the low-level numeric routines, you will be fighting numeric garbage as some of the values involved are easily more than 2^53 times other values and thus proper cancellation cannot occur when using iterative floating point methods.

请先登录，再进行评论。

类别

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

标签

Community Treasure Hunt

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

Start Hunting!

Translated by