Info

Questa domanda è chiusa. Riaprila per modificarla o per rispondere.

why a variable is remaining unchanged when solving two simultaneous ordinary differential equation using ode45?

1 visualizzazione (ultimi 30 giorni)
function dX = example1(t,x)
To=400;Tc=350; H = -47.8012;
k=1.97*10^20; E=166*10^3; V = 3000;
Ca0=9; m=4; rho= 1; UA=7;
Cp=1;R=8.314;
dX=[(m*(Ca0-x(1))/rho*V)-k*x(1)*exp(-E/R*x(2));
((m*Cp*(To-x(2))-(V*H*x(1)*k*exp(-E/R*x(2)))+UA*(Tc-x(2)))/V*rho*Cp)];
in Command window:
>>tspan = [0 1000];x0=[9 400];
>>[t,x] = ode45(@example1,tspan,x0);
In the output of x the value of x(1) remained unchanged throughout the timespan where the x(2) was changed. I tried changing time span to 10^6 and furthur higher. This did not solve the problem. Can anyone kindly suggest the mistake i had done or changes to the program so that x(1) also changes simultaneously.

Risposte (1)

Steven Lord
Steven Lord il 29 Mag 2020
Let's look at both terms of the expression you use to define dX(1).
dX=[(m*(Ca0-x(1))/rho*V)-k*x(1)*exp(-E/R*x(2));
For the first term, (m*(Ca0-x(1))/rho*V), the initial value of x(1) is 9 and so is Ca0. This means this term is 0 at the initial conditions.
In the second term, the quantity E/R is approximately 20000. When you compute exp(-E/R*x(2)), unless x(2) is less than about 1/30 this quantity will underflow to 0. At the initial conditions, x(2) is 400.
You may say "But wait, k is pretty big! Doesn't that counteract the small value of that exponential?" If we were to compute this symbolically to see just how small that part of this term is:
>> vpa(exp(-E/R*sym(400)))
ans =
4.3226475784249548860556479944784e-3468506
Yeah, that's zero in double precision. Even if you were to multiply it by k in arbitrary precision before converting back to double precision, it's still zero.
So at your initial conditions, dX(1) is 0. This means x(1) doesn't change from its initial value. At the next time step, the first term in that expression remains 0. I haven't examined the part of the expression for dX that defines dX(2), but I'm guessing x(2) doesn't decrease down below 1/30 all that quickly so that term probably also remains 0. So at the next time step, x(1) probably remains the same.
Are you certain you're using the correct equations and the correct units for the parameters? Try doing dimensional analysis to check you're not writing E in units of millimeters and R in units of meters or something similar.
  4 Commenti
Steven Lord
Steven Lord il 29 Mag 2020
You're probably going to need to ask your professor or teaching assistant for help understanding why your code isn't getting the same answer as your textbook. We could blindly start tweaking the values of your constants or your equations, but doing so isn't likely to satisfy the people grading your assignment.

Prodotti


Release

R2016a

Community Treasure Hunt

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

Start Hunting!

Translated by