3 views (last 30 days)

Many greetings to anyone who stumbles upon this thread. I'm currently trying to solve this differential equation for a chemical reaction rate regression analysis.

As you can see, I have 15 parameters that I've already fitted and real & known values. The differential equation has DCF as y and O3 and FeOOH are also known constant values. My only problem is writing a code that continously finds the roots of [.OH] to input inside the differential equation. My time span is 20 minutes and my initial value of DCF is 32.439. I've tried with the following code:

function dy = ozecat(t,y,x,g,ox)

OH = [x(4) (x(5)*y+x(6)*g(1)+x(7)*g(1)/y-x(8)*ox(1)) (x(9)*ox(1)*y+x(10)*(ox(1)^2)-x(11)*ox(1)*g(1)+x(12)*(ox(1)*g(1)/y)) (-x(13)*ox(1)^2*y-x(14)*ox(1)^2*g(1)+x(15)*(ox(1)^2*g(1)/y))];

dy = -x(1)*ox(1)*y-x(2)*y*roots(OH)-x(3)*g(1)*ox(1);

end

where g stands for FeOOH data, ox for O3 data and x are my array of 16 parameters (from a to l and the k's). I've also tried to select one single root of OH by adding this two lines of code to the above:

rroot_OH = roots(OH);

rroot_OH = rroot_OH(imag(rroot_OH) == 0;

But to no avail. I get the following errors. This first for the original code:

Error using horzcat

Dimensions of matrices being concatenated are not consistent.

Error in ozecat (line 2)

OH = [x(4) (x(5)*y+x(6)*g(1)+x(7)*g(1)/y-x(8)*ox(1)) (x(9)*ox(1)*y+x(10)*(ox(1)^2)-x(11)*ox(1)*g(1)+x(12)*(ox(1)*g(1)/y))

(-x(13)*ox(1)^2*y-x(14)*ox(1)^2*g(1)+x(15)*(ox(1)^2*g(1)/y))];

Error in odearguments (line 90)

f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.

and this one for the modified with the two lines added:

Error using odearguments (line 95)

OZECAT returns a vector of length 3, but the length of initial conditions vector is 1. The vector returned by OZECAT and the initial

conditions vector must have the same number of elements.

Error in ode45 (line 115)

odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in parametrows (line 34)

[t_calc y] = ode45(@ozecat, tspan, y0,[],x,g,ox);

So how exactly is this done? Or is it not possible by ode45?

Steven Lord
on 4 Dec 2019

This doesn't look like a system of ODEs. The presence of the variable d in the first set of curly braces threw me for a second, but I think this is a system of one differential equation and one algebraic equation. That makes it a system of differential algebraic equations (DAEs.) You may want to use the technique shown for the Robertson problem on that page as a model for your DAE function.

If you do this, you may want to choose names for your variables inside your DAE function carefully so they look closer to your mathematical equations. This may help you if you need to debug or present your work to someone more familiar with chemical reaction rate regression analysis than I am.

function dDCF = ozecat(t,DCF,params,FeOOH,O3)

a = params(1); % params used to be called x

b = params(2);

% etc

% Not sure how OH is calculated

OH = ...;

term(1) = a*OH.^3;

term(2) = (OH.^2).* ...;

% ...

dDCF = ...;

end

Sign in to answer this question.

Opportunities for recent engineering grads.

Apply Today
## 3 Comments

## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/494867-how-to-solve-a-ode-that-is-a-function-of-a-cubic-polynomial-f-x-0#comment_774501

⋮## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/494867-how-to-solve-a-ode-that-is-a-function-of-a-cubic-polynomial-f-x-0#comment_774501

## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/494867-how-to-solve-a-ode-that-is-a-function-of-a-cubic-polynomial-f-x-0#comment_774506

⋮## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/494867-how-to-solve-a-ode-that-is-a-function-of-a-cubic-polynomial-f-x-0#comment_774506

## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/494867-how-to-solve-a-ode-that-is-a-function-of-a-cubic-polynomial-f-x-0#comment_774520

⋮## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/494867-how-to-solve-a-ode-that-is-a-function-of-a-cubic-polynomial-f-x-0#comment_774520

Sign in to comment.