ODE and Data fitting

I have 3 differential equations:

The initial condition for all the equations at timepoint=0 equal to 0. After solving this equaitons, I want to a data fitting on the equation below in order to find the nest fitted A, B, and C.

Can someone help me on this, because I face several errors when I tried to do it.

Maybe the best way is to find the u, v, and w functions and then place them in equations x and y; but how can I do it? is there any way that ode45 give the formula of u,v, and w?

Answer 1

[Star Strider]

All I can do is to point you in the direction of successful efforts to do what you want.

One such is: Parameter Estimation for a System of Differential Equations and another is of course Monod kinetics and curve fitting.

We can’t help without seeing your code.

[Ehsan Homaee]

Thanks for your reply. I atached my code to this. I think my problem differs from you mentioned. Could you please guide me what should I do in this case?

[Star Strider]

First, do not use global variables! They create code that is almost impossible to debug.

Second, you apparently have three differential equations (I am going on your original Question, since I cannot decipher your code) and you are fitting one variable. I have no idea which one that is.

My ‘Monod kinetics’ code does something similar, so you can use that as a prototype.

[Star Strider]

The updated post clarified everything. This runs and gives a reasonably decent fit, however it may be necessary to use the optimoptions function to provide an options structure to allow lsqcurvefit a higher function evaluation limit. I assume the parameters are allowed to become negative (not always permitted), so I did not constrain them, since you did not in your code.

Try this:

    function Crodata1_Fit
        function y = Objfcn(b,t)
            x0 = b(end-2:end);
            
            [tv,xv] = ode45(@ODEs, t, x0);
            
            function dY = ODEs(T,Y)                     % Y(1) = v;  Y(2) = u,  Y(3) = w
                dY(1,:) = Y(1).*(-1.0./2.0e+1)+1.0./(Y(1).^2.*2.5e+1+Y(2).^2.*2.5e+1+2.5e+1);
                dY(2,:) = (Y(3).^2.*(3.0./5.0e+1))./(Y(3).^2+1.0)+(Y(2).^2.*(2.0./2.5e+1))./(Y(1).^2+Y(2).^2+1.0)+Y(2)./2.5e+1;
                dY(3,:) = Y(3).*(-3.0./2.5e+1)+1.0./(Y(1).^2.*2.5e+1+Y(2).^2.*2.5e+1+2.5e+1);
            end
            y = 1./(b(1) + b(2).*xv(:,2).^2 + b(3).*xv(:,1).^2);
        end
        
        tdata=[0	0.5	1	2	3	5	7.5	10	20	30	60	120	180].';
        Crodata1=[0.037373477
            2.938766958
            2.428716866
            4.673151154
            4.039997314
            7.286068277
            5.229206437
            6.722726692
            1.809539197
            1.975990902
            1.366188561
            1.103164147
            0.29209409];
        
        B0 = [rand(3,1); zeros(3,1)];
        
        [B,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@Objfcn,B0,tdata,Crodata1);
        
        fprintf(1,'\tRate Constants:\n')
        for k1 = 1:length(B)
            fprintf(1, '\t\tB(%d) = %8.5f\n', k1, B(k1))
        end
        
        tv = linspace(min(tdata), max(tdata));
        Yfit = Objfcn(B, tv);
        
        figure(1)
        plot(tdata, Crodata1, 'p')
        hold on
        hlp = plot(tv, Yfit);
        hold off
        grid
        xlabel('Time')
        ylabel('Unknown Unit')
    end

I run all this code inside of a specific function I use for such purposes, so to make life easier, I named this one ‘Crodata1_Fit’. Save it as a separate function and run it by calling:

    Crodata1_Fit

from a script or your Command Window.

Make appropriate changes to get the result you want.

One run produced this plot:

and these parameters:

		B(1) =  0.14574
		B(2) =  0.00926
		B(3) =  2.80074
		B(4) = -0.32271
		B(5) = -0.02885
		B(6) = -1.57333

Other optimisers (specifically the ga funciton) will likely produce even better results.

NOTE —

I used this code to derive the ‘ODEs’ equations:

syms u(t) v(t) w(t) t T Y 
DEs = [diff(u)==0.08*u^2/(1+u^2+v^2)+0.06*w^2/(1+w^2)+0.04*u; diff(v)==0.04*1/(1+u^2+v^2)-0.05*v; diff(w)==0.04*1/(1+u^2+v^2)-0.12*w];
[VF,Sbs] = odeToVectorField(DEs)
ODEs = matlabFunction(VF, 'Vars',{T,Y})

It is easier (and more accurate) to let the Symbolic Math Toolbox do the ‘heavy lifting’.

[Star Strider]

Ehsan Homaee’s Answer moved here —

Thanks for your response. Why the output has 6 diffeent B? Since the objective function only has 3 parameters,what is the other 3 parameters?

Furthermore, if I want to restrict the parameters to only possitive values, how can I do it?

[Star Strider]

As always, my pleasure!

The first 3 ‘B’ parameters correspond to ‘A’, ‘B’, and ‘C’. The last 3 are the initial conditions of the system of differential equations, since I always let the optimization function estimate them as well.

To restrict the first 3 to positive values, the lsqcurvefit call changes to:

[B,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@Objfcn,B0,tdata,Crodata1, [zeros(1,3), -Inf(1,3)]);

I would not constrain the initial conditions. However if you want to constrain them as well to be positive:

[B,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@Objfcn,B0,tdata,Crodata1,zeros(1,6));

I suspect the ga function would come up with the best parameter set, since it searches the entire parameter space and is not usually affected by local minima.