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

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

[Walter Roberson]

Correspondance: initial conditions are:

A(0) = 10

B(0) = 100

C(0) = 0

Answer 1

[Matt Tearle]

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.

Answer 2

[James]

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

[Matt Tearle]

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.

Answer 3

[Walter Roberson]

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))

[Walter Roberson]

I thought I had a symbolic solution, but I had a typo early on :(

[Walter Roberson]

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]

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.