How to solve a ODE that is a function of a cubic polynomial f(x) = 0?
1 visualizzazione (ultimi 30 giorni)
Mostra commenti meno recenti
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?
3 Commenti
Risposte (1)
Steven Lord
il 4 Dic 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
Vedere anche
Categorie
Scopri di più su Ordinary Differential Equations in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!