Can the equations in the differential equation group be represented by symbols when matlab ode solves the differential equation group ?

main function:

function dcdzh2h2sco2_centri = dcdzh2h2sco2_centri(z,c)
R = 8.314;
g = 1000;
Mi_h2s = 34;
Mi_co2 = 44;
Mi_ch4 = 16;
syms T z c_ch4 c_co2 c_h2s dc_31 dc_32 dc_33
Tz = 293 + 10*z;
alpha_1_2 = log((3.6029/(1.3932))*((c_co2)/c_h2s))/(log(Tz/293.15));
alpha_1_3 = log((3.6029/(37.063))*((c_ch4)/c_h2s))/(log(Tz/293.15));
alpha_2_3 = log((1.3932/(37.063))*((c_ch4)/c_co2))/(log(Tz/293.15));
alpha_h2sT = (c_co2/(c_ch4+c_co2+c_h2s))*(-alpha_1_2)+((c_ch4)/(c_ch4+c_co2+c_h2s))*(-alpha_1_3);
alpha_co2T = (c_ch4/(c_ch4+c_co2+c_h2s))*(-alpha_2_3)+((c_h2s)/(c_ch4+c_co2+c_h2s))*(alpha_1_2);
alpha_ch4T = (c_h2s/(c_ch4+c_co2+c_h2s))*(alpha_1_3)+((c_co2)/(c_ch4+c_co2+c_h2s))*(alpha_2_3);
alpha_h2s = subs(alpha_h2sT,T,Tz);
alpha_co2 = subs(alpha_co2T,T,Tz);
alpha_ch4 = subs(alpha_ch4T,T,Tz);
Hi_h2s = 1/((Mi_h2s*g)/(R*Tz));
Hi_co2 = 1/((Mi_co2*g)/(R*Tz));
Hi_ch4 = 1/((Mi_ch4*g)/(R*Tz));
dcdz_h2s_1 = -(0.027*((1 + alpha_h2s)/Tz)-(1/Hi_h2s))
dcdz_co2_1 = (-(0.027*((1 + alpha_co2)/Tz)-(1/Hi_co2)))
dcdz_ch4_1 = (-(0.027*((1 + alpha_ch4)/Tz)-(1/Hi_ch4)))
dc_11 = char(dcdz_h2s_1)
dc_12 = char(dcdz_co2_1)
dc_13 = char(dcdz_ch4_1)
dc_21 = strrep(dc_11,'c_h2s','c(1)');
dc_21 = strrep(dc_21,'c_co2','c(2)');
dc_21 = strrep(dc_21,'c_ch4','c(3)')
dc_22 = strrep(dc_12,'c_h2s','c(1)');
dc_22 = strrep(dc_22,'c_co2','c(2)');
dc_22 = strrep(dc_22,'c_ch4','c(3)')
dc_23 = strrep(dc_13,'c_h2s','c(1)');
dc_23 = strrep(dc_23,'c_co2','c(2)');
dc_23 = strrep(dc_23,'c_ch4','c(3)')
dc_31 = str2sym(dc_21)
dc_32 = str2sym(dc_22)
dc_33 = str2sym(dc_23)
dcdzh2h2sco2_centri = zeros(3,1);
dcdzh2h2sco2_centri = [dc_31;dc_32;dc_33]
end

and solve it:

[z,c] = ode45(@dcdzh2h2sco2_centri,[0:0.1:100],[3.6029 1.3932 37.063])

it gives results and errors as following:

dcdzh2h2sco2_centri =

((27*c(3)*log((7004729934415739*c(3))/(72057594037927936*c(1))))/(log((200*z)/5863 + 5860/5863)*(c(1) + c(2) + c(3))) + (27*c(2)*log((36029*c(2))/(13932*c(1))))/(log((200*z)/5863 + 5860/5863)*(c(1) + c(2) + c(3))) - 27)/(10000*z + 293000) + 1/((4157*z)/1700000 + 1218001/17000000)

1/((4157*z)/2200000 + 1218001/22000000) - ((27*c(1)*log((36029*c(2))/(13932*c(1))))/(log((200*z)/5863 + 5860/5863)*(c(1) + c(2) + c(3))) - (27*c(3)*log((5417297035514729*c(3))/(144115188075855872*c(2))))/(log((200*z)/5863 + 5860/5863)*(c(1) + c(2) + c(3))) + 27)/(10000*z + 293000)

1/((4157*z)/800000 + 1218001/8000000) - ((27*c(2)*log((5417297035514729*c(3))/(144115188075855872*c(2))))/(log((200*z)/5863 + 5860/5863)*(c(1) + c(2) + c(3))) + (27*c(1)*log((7004729934415739*c(3))/(72057594037927936*c(1))))/(log((200*z)/5863 + 5860/5863)*(c(1) + c(2) + c(3))) + 27)/(10000*z + 293000)

Error using odearguments ( line 113 )

Input must be single precision or double precision floating point value.

I want to use the symbols deducing results to represent diffenrential equation in the differential equation groups,could you help me?Thanks a lot.

Answer 1

[Torsten]

function dcdzh2h2sco2_centri = dcdzh2h2sco2_centri(z,c)  
R = 8.314;
g = 1000;
Mi_h2s = 34;
Mi_co2 = 44;
Mi_ch4 = 16;
syms T z c_ch4 c_co2 c_h2s dc_31 dc_32 dc_33 Z
Tz = 293 + 10*z;
alpha_1_2 = log((3.6029/(1.3932))*((c_co2)/c_h2s))/(log(Tz/293.15));
alpha_1_3 = log((3.6029/(37.063))*((c_ch4)/c_h2s))/(log(Tz/293.15));
alpha_2_3 = log((1.3932/(37.063))*((c_ch4)/c_co2))/(log(Tz/293.15));
alpha_h2sT = (c_co2/(c_ch4+c_co2+c_h2s))*(-alpha_1_2)+((c_ch4)/(c_ch4+c_co2+c_h2s))*(-alpha_1_3);
alpha_co2T = (c_ch4/(c_ch4+c_co2+c_h2s))*(-alpha_2_3)+((c_h2s)/(c_ch4+c_co2+c_h2s))*(alpha_1_2);
alpha_ch4T = (c_h2s/(c_ch4+c_co2+c_h2s))*(alpha_1_3)+((c_co2)/(c_ch4+c_co2+c_h2s))*(alpha_2_3);
alpha_h2s = subs(alpha_h2sT,T,Tz);
alpha_co2 = subs(alpha_co2T,T,Tz);
alpha_ch4 = subs(alpha_ch4T,T,Tz);
Hi_h2s = 1/((Mi_h2s*g)/(R*Tz));
Hi_co2 = 1/((Mi_co2*g)/(R*Tz));
Hi_ch4 = 1/((Mi_ch4*g)/(R*Tz));
dcdz_h2s_1 = -(0.027*((1 + alpha_h2s)/Tz)-(1/Hi_h2s));
dcdz_co2_1 = (-(0.027*((1 + alpha_co2)/Tz)-(1/Hi_co2)));
dcdz_ch4_1 = (-(0.027*((1 + alpha_ch4)/Tz)-(1/Hi_ch4)));
f = matlabFunction([dcdz_h2s_1;dcdz_co2_1;dcdz_ch4_1]);
dcdzh2h2sco2_centri = f(c(3),c(2),c(1),z);
end

Doublecheck whether the order of the arguments for f is correct, i.e. look if f has the form

@(c_ch4, c_co2, c_h2s, z) [f1;f2;f3]

(I guess c_ch4 corresponds to c(3), c_co2 corresponds to c(2) and c_h2s corresponds to c(1) in your code)

[Hou X.Y]

Thans a lot for your answering,but it still throws errors:

Error using odearguments ( line 113 )

Input must be single precision or double precision floating point value.

when I use

[z,c] = ode45(@dcdzh2h2sco2_centri,[0:0.1:100],[3.6029 1.3932 37.063])

solve it.

[Torsten]

How does f look like after the command "f = matlabFunction(...)" ?

Note also that you do this conversion each time "dcdzh2h2sco2_centri" is called. This costs a lot of time.

You should generate the function handle in the program from where you call ode45 and already supply it to ode45 in its call:

% Your calculations
f = matlabFunction(...);
F = @(z,c) f(c(3),c(2),c(1),z);
[T,Y] = ode45(F,...)

[Hou X.Y]

Thank you very much.It does work!

[Torsten]

My pleasure.

Answer 2

[John D'Errico]

Sorry, but the numerical ode solvers, such as ODE45 are NOT symbolic tools. They CANNOT solve an ode with symbolic parameters. Period.

If you want to solve an ode using symbolic parameters then you would need to use a tool like DSOLVE. Of course, there is no assurance a symbolic solution can be found. It would be my guess that most such complicated odes will fail in that respect, lacking any symbolic solution, but you may always get lucky.

Such is life in the big city. There are many things in life we may want, but cannot have.

[Hou X.Y]

Thanks for your answer.But it seems that you have a wrong understanding of my question:what I want is not use ode to get symbols solution,I want to let the symbols(eg.dc_31) represent the results of symbolic operation,then let them represent the equations in a differential equation group,then use ode to solve it and get numerical solution.