Ode system solution with unknown constant

I'm trying to obtain and plot a function that adjusts to my data. I have a table with time, position and velocity, and an ode function which I have rewritten as a system:

ydot(1) = y(2)

ydot(2) = -a*y(2) - b*y(1) + c*cos(wt)

a,b,c and w are unknown parameters but I can estimate their aproximate values.

How can I solve the problem?

Answer 1

[Star Strider]

You have already seen and posted a Comment to Parameter Estimation for a System of Differential Equations and that (or similar threads) are what I would direct you to.

For your application, the ‘kinetics’ objective function in that solution would be:

function Y = kinetics(theta,t)
a = theta(1);
b = theta(2);
c = theta(3);
w = theta(4);
y0 = theta(5:6);

[t,yv] = ode45(@Difeq,t,y0);

function ydot = Difeq(t,y)
ydot = zeros(2,1);
ydot(1) = y(2);
ydot(2) = -a*y(2) - b*y(1) + c*cos(w*t);
end

Y = yv;
end

Note that here the code will estimate the intial conditions as well, so ‘theta0’ must be a 6-element vector.

I obviously cannot test this, so I am posting it as UNTESTED CODE. It should work, although it may require a bit of editing.

[ana ramirez de las heras]

Thank you so much!! I will test it and try to make it work! You are very helpful, I will let you know if I achieve the solutions!

[Star Strider]

My pleasure!

If I had your data, I could test it.

[ana ramirez de las heras]

Hello Star Strider!

The code is working! Find the data (datos.txt) and the .m file attached. I found out that the a,b,c and w values that I had estimated were not so accurate...

The columns correspond to time, position and velocity. The ode function which represents the data is

y'' + ay' + by = c cos (wt)

Thank you for your help!!

[Star Strider]

My pleasure!

I could not get a good fit with lsqcurvefit, so I went with ‘Plan B’ and wrote ga (genetic algorithm) code for it. I added a phase term to the cos call (entirely appropriate, and does not assume that the cos term always begins at 0) and while the fit improved slightly, it did so only for position (matching the frequency and phase, however not the amplitude), not velocity, that remains difficult (essentially impossible) to fit.

Using this ‘solver’ code with the added phase term:

%Calcular los parametros que se ajustan a los datos
function Y=solver(theta,t)
a = theta(1);
b = theta(2);
c = theta(3);
w = theta(4);
p = theta(5);                                       % Phase
y0 = theta(6:7);
[t,yv] = ode45(@Difeq,t,y0);
    function ydot = Difeq(t,y)
    ydot = zeros(2,1);
    ydot(1) = y(2);
    ydot(2) = -a*y(2) - b*y(1) + c*cos(w*t + p);
    end
Y = yv;
end

likely the best parameter set ga could estimate were:

	Rate Constants:
		Theta(1) =  4.22996
		Theta(2) =  6.33071
		Theta(3) = 19.42851
		Theta(4) =  2.01714
		Theta(5) =  5.17042
		Theta(6) =  0.69420
		Theta(7) =  2.37232

with a fitness value of 63.9.

If you have the Global Optimization Toolbox, I will post my ga code for your problem.

[Star Strider]

My pleasure!

I discovered that the columns of ‘y’ are incorrect, and should be reversed, the correct orientation being:

y=[v p]; 

since the derivatrive is the first column of ‘Y’ and its integral is the second column. (I was too tired yesterday to detect that error in your code.)

Make that change, and the fit is excellent:

	Rate Constants:
		Theta(1) =  0.04702
		Theta(2) =  9.99643
		Theta(3) = 20.76383
		Theta(4) =  1.99939
		Theta(5) =  3.17972
		Theta(6) =  0.94858
		Theta(7) =  1.39551

The fitness value (norm of the residuals) for this run was typical at 13.6113.

I plotted the two variables separately, since it is difficult to see them when all are plotted together.

[Star Strider]

As always, my pleasure!

I typically let the optimisation function calculate everything, in order not to constrain it beyond any obvious limits (for example constraining ‘w’ and ‘p’ to here).