Error on dsolve() line 201. Solving system of differential equations.

Long post warning.

Hi there. I'm new on Matlab. I've been watching YouTube videos for a week and I can't figure out the right way to solve an engineering problem. It's a Engeneering of chemical reactions problem.

Here are some screenshots of the problem, it was supposedly already solved on Polymath.

Some initial conditions here.

Here is the solution.

And here's what my teacher did on Polymath.

I want to learn how to use this tool properly and that's why I'm asking this. I have the following code I wrote myself, but it doesn't work.

%Se debe limpiar la pantalla
clear all;
clc;

%Se definen las variables que dependen de otras
%Commented out this because I'm not sure if they are needed
%syms kA(T) CA(x,y,T) CB(x,y,T) CC(X,y,T) kC(T) x(W) y(W) T(W);
syms x(W) y(W) T(W);

%Definir variables que no cambian... o sea, constantes
T0=325;
TA=325;
Ea=25000;
FA0=5;
FB0=FA0;
FI=FA0;
k=0.004; %a 310 K
Ke=1000; %a 303 K
DHRx=-20000;
alfa=0.0002;
Cp=20; CA=Cp; CB=Cp; CC=Cp; CpI=40;
UaRoB=0.5;
R=1.987;
Xeq=(Ke^(1/2))/(2+Ke^(1/2)); %Se obtiene con -rA = 0
y0=1/3; Ctot=0.3; CA0=y0*Ctot;

%Se declaran variables que dependen de otras
kA=k*exp((Ea/R)*(1/310-1/T));
kC=Ke*exp((DHRx/R)*(1/303-1/T));
CA=CA0*(1-x)*y*(T0/T);
CB=CA0*(1-x)*y*(T0/T);
CC=2*CA0*x*y*(T0/T);
rA=-kA*(CA*CB-((CC^2)/kC));
sCp=CA+CB+CC;

%Definir ecuaciones diferenciales
ode1 = diff(x) == -1*rA/FA0;
ode2 = diff(y) == -(alfa/2*y)*(T/T0);
ode3 = diff(T) == (UaRoB*(TA-T)+(-DHRx*-rA))/(FA0*sCp);
odes = [ode1;ode2;ode3];

%Se almacena solucion en una variable
Solucion = dsolve(odes);

%Extrae la función del arreglo de la solución
Solx(W) = Solucion.x
Soly(W) = Solucion.y
SolT(W) = Solucion.T

[Solx(W), Soly(W), SolT(W)] = dsolve(odes);

cond1 = x(0) == 0;
cond2 = y(0) == 1;
cond3 = T(0) == 325;
conds = [cond1; cond2; cond3];

[Solx(W), Soly(W),SolT(W)] = dsolve(odes,conds);
plot(y,W(:,1))
xlabel('Peso de catalizador (W)'),ylabel('y')

There is a problem on the solver, dsolve at line 201. On a help page from this site, I stumbled upong a very clear and easy-to-follow routine, but it doesn't seem to work with my specific case.

How can I improve the code and get the solutions? I feel dumb, but it's my first time dealing with this software.

Answer 1

[Torsten]

Try

function main
z0=[0; 1; 325];
tspan = [0 2500];
[T Z]=ode15s(@myfun,tspan,z0);
plot(T,Z(:,1),'-o',T,Z(:,2),'-.',T,Z(:,3),'-*')

function dz = myfun(W,z)  
x=z(1);
y=z(2);
T=z(3);

dz=zeros(3,1);

%Definir variables que no cambian... o sea, constantes
T0=325;
TA=325;
Ea=25000;
FA0=5;
FB0=FA0;
FI=FA0;
k=0.004; %a 310 K
Ke=1000; %a 303 K
DHRx=-20000;
alfa=0.0002;
Cp=20; CA=Cp; CB=Cp; CC=Cp; CpI=40;
UaRoB=0.5;
R=1.987;
Xeq=(Ke^(1/2))/(2+Ke^(1/2)); %Se obtiene con -rA = 0
y0=1/3; Ctot=0.3; CA0=y0*Ctot;

%Se declaran variables que dependen de otras
kA=k*exp((Ea/R)*(1/310-1/T));
kC=Ke*exp((DHRx/R)*(1/303-1/T));
CA=CA0*(1-x)*y*(T0/T);
CB=CA0*(1-x)*y*(T0/T);
CC=2*CA0*x*y*(T0/T);
rA=-kA*(CA*CB-((CC^2)/kC));
sCp=CA+CB+CC;

%Definir ecuaciones diferenciales
dz(1) = -1*rA/FA0;
dz(2) = -(alfa/2*y)*(T/T0);
dz(3) = (UaRoB*(TA-T)+(-DHRx*-rA))/(FA0*sCp);

Best wishes

Torsten.

[Joseluis Colunga]

It worked! You are a genius!

Care to explain why did you change the structure that drastically?

[Torsten]

I don't think that the structure has changed that much.

I only replaced symbolic to the usual numeric structures.

Best wishes

Torsten.

Answer 2

[Joseluis Colunga]

Anyone?

[Torsten]

You won't succeed using "dsolve" because your system of ODEs is too complicated.

Use a numerical solver instead, e.g. ODE15S.

Best wishes

Torsten.

[Joseluis Colunga]

Thank you very much for your help.

Answer 3

[Walter Roberson]

The information I find about PolyMath 5.1 suggests that it is a numeric differential equations solver.

MATLAB's dsolve() is an analytic solver. It appears that dsolve() is quickly giving up trying to find an analytic solution to the equations, and so is returning an empty symbolic variable.

Your second equation forces y(W) to be exponential in form, and your other equations then force it those terms to be at least exponential themselves. The system is messy. I am not sure there is an analytic solution.

You should consider using one of the numeric solvers such as ode45().

I seem to remember that other people working with chemical reactions have had to switch to stiff solvers such as ode15s.

[Joseluis Colunga]

But can I solve them simultaneously? Storing the diff eq. on an array? Like the one from:

odes = [ode1;ode2;ode3];

I have to get these graphs:

And I realized I had to assign values to W with

W=[0 2500];

Since it is the size of the x axis from my graphs, quite an arbitrary number, but W is already defined.

At the end, it all ends up depending on x, y, T. Or every variable on each other.

[Walter Roberson]

Create one function, with time as the first input (I think), and a vector of values as the second input. That vector of values should probably be 6 elements long: one each for x, y, and T, and one each for dx, dy, and dT .

.... something like that. I am a little fuzzy on exactly what inputs to use for the function.

[Joseluis Colunga]

Thank you very much for your help.