Help with solving systems differential equations

Hello everyone,

I have the following set of coupled first-order differential equations:

dA/A+dB/B=0

-dC-dD/kn1=R*A*dA

kn2*dE=-A*dA

dC/C=dE/E+dB/B

ds=kn2*dE/E-kn3*dC/C

dh=kn2*dE

dD/kn1=kn4*B*A^2/2/kn5

Where kn1,....kn5 are known parametres

I tried using the method that use here https://es.mathworks.com/matlabcentral/answers/514307-solve-numerically-a-system-of-first-order-differential-equations but matlab say it's unable to result it

I was wondering which could be a good attempt to solve numerically this system of differential equations.

Any suggestion?

Answer 1

[Walter Roberson]

syms A(t) B(t) C(t) D(t) E(t) h(t) s(t)
syms kn1 kn2 kn3 kn4 kn5 R
dA = diff(A);
dB = diff(B);
dC = diff(C);
dD = diff(D);
dE = diff(E);
ds = diff(s);
dh = diff(h);

eqn1 = dA/A + dB/B == 0

eqn1(t) = 

eqn2 = -dC - dD/kn1 == R * A * dA

eqn2(t) = 

eqn3 = kn2*dE == -A*dA 

eqn3(t) = 

eqn4 = dC/C == dE/E + dB/B

eqn4(t) = 

eqn5 = ds == kn2 * dE/E - kn3 * dC/C

eqn5(t) = 

eqn6 = dh == kn2*dE

eqn6(t) = 

eqn7 = dD/kn1 == kn4*B*A^2/2/kn5

eqn7(t) = 

eqns = [eqn1,eqn2,eqn3,eqn4,eqn5,eqn6,eqn7].'

eqns(t) = 

fcns = [A(t),B(t),C(t),D(t),E(t),h(t),s(t)];
[M, F] = massMatrixForm(eqns, fcns)

M = 

F = 

%convert matrix form into function of variables
syms t x [1 7]
Q = subs(M*ones(7,1)-F, fcns, x)

Q = 

%parameter values for demonstration purposes
kn_vals = rand(1,5)

kn_vals = 1×5

0.8901 0.1594 0.4955 0.6421 0.5689

R_vals = rand*10

R_vals = 7.3070

%create a function handle
obj = matlabFunction(subs(Q,[kn1,kn2,kn3,kn4,kn5,R], [kn_vals,R_vals]), 'vars', {t, x(:)})

obj = function_handle with value:
    @(in1,in2)[1.0./in2(1,:)+1.0./in2(2,:);in2(1,:).*(-7.306962374164408)-2.123443587802794;in2(1,:)+1.594441238491097e-1;-1.0./in2(2,:)+1.0./in2(3,:)-1.0./in2(5,:);4.95537979655735e-1./in2(3,:)-1.594441238491097e-1./in2(5,:)+1.0;8.405558761508903e-1;in2(1,:).^2.*in2(2,:).*(-5.643850562679325e-1)+1.123443587802794]

%initial conditions for demonstration purposes
x0 = rand(1,7)

x0 = 1×7

0.7536 0.2815 0.8073 0.9225 0.3860 0.4098 0.0463

%run the function
[t,X] = ode45(obj, [0 0.015], x0);
plot(t,X)

For the random coefficients that I put in, you cannot get much further than this; shortly after this point, it fails integration, and that happens for both ode45() and ode23s()

[SARA TORRES UGIA]

Hi, thanks for your answer. I try do it with the real values and got this error:

Error using sym/subs>normalize (line 240)
Entries in second argument must be scalar.

Error in sym/subs>mupadsubs (line 166)
[X2,Y2,symX,symY] = normalize(X,Y); %#ok

Error in sym/subs (line 154)
    G = mupadsubs(F,X,Y);

Error in resolver (line 40)
obj = matlabFunction(subs(Q,[kn1,kn2,kn3,kn4,kn5], [kn_vals]), 'vars', {t, x(:)})

This is the code that I use:

syms A(t) B(t) C(t) D(t) E(t) h(t) s(t)
syms kn1 kn2 kn3 kn4 kn5 R
dA = diff(A);
dB = diff(B);
dC = diff(C);
dD = diff(D);
dE = diff(E);
ds = diff(s);
dh = diff(h);

eqn1 = dA/A + dB/B == 0
eqn2 = -dC - dD/kn1 == B * A * dA
eqn3 = kn2*dE == -A*dA 
eqn4 = dC/C == dE/E + dB/B
eqn5 = ds == kn2 * dE/E - kn3 * dC/C
eqn6 = dh == kn2*dE
eqn7 = dD/kn1 == kn4*B*A^2/2/kn5

eqns = [eqn1,eqn2,eqn3,eqn4,eqn5,eqn6,eqn7].'

fcns = [A(t),B(t),C(t),D(t),E(t),h(t),s(t)];
[M, F] = massMatrixForm(eqns, fcns)

%convert matrix form into function of variables
syms t x [1 7]
Q = subs(M*ones(7,1)-F, fcns, x)

%parameter values for demonstration purposes
kn1= 0.001963;
kn2= 1005;
kn3= 287;
kn4= 0.023;
kn5= 0.05;
kn_vals = [kn1 kn2 kn3 kn4 kn5]

%create a function handle
obj = matlabFunction(subs(Q,[kn1,kn2,kn3,kn4,kn5], [kn_vals]), 'vars', {t, x(:)})

%initial conditions for demonstration purposes
x0 = [85 1.703 220000 0 450 0 0];

%run the function
[t,X] = ode45(obj, [0 0.0005], x0);
plot(X,t)

[Walter Roberson]

syms A(t) B(t) C(t) D(t) E(t) h(t) s(t)
syms kn1 kn2 kn3 kn4 kn5 R
dA = diff(A);
dB = diff(B);
dC = diff(C);
dD = diff(D);
dE = diff(E);
ds = diff(s);
dh = diff(h);

eqn1 = dA/A + dB/B == 0
eqn2 = -dC - dD/kn1 == B * A * dA
eqn3 = kn2*dE == -A*dA 
eqn4 = dC/C == dE/E + dB/B
eqn5 = ds == kn2 * dE/E - kn3 * dC/C
eqn6 = dh == kn2*dE
eqn7 = dD/kn1 == kn4*B*A^2/2/kn5

eqns = [eqn1,eqn2,eqn3,eqn4,eqn5,eqn6,eqn7].'

fcns = [A(t),B(t),C(t),D(t),E(t),h(t),s(t)];
[M, F] = massMatrixForm(eqns, fcns)

%convert matrix form into function of variables
syms t x [1 7]
Q = subs(M*ones(7,1)-F, fcns, x)

%parameter values for demonstration purposes
Kn1= 0.001963;
Kn2= 1005;
Kn3= 287;
Kn4= 0.023;
Kn5= 0.05;
kn_vals = [Kn1 Kn2 Kn3 Kn4 Kn5]

%create a function handle
obj = matlabFunction(subs(Q,[kn1,kn2,kn3,kn4,kn5], [kn_vals]), 'vars', {t, x(:)})

%initial conditions for demonstration purposes
x0 = [85 1.703 220000 0 450 0 0];

%run the function
[t,X] = ode45(obj, [0 0.0005], x0);
plot(t,X)