Numerically solve a 3rd order differential equation with 3 unknowns

Hi, I have a 3rd order differential equation in y, with 3 unknowns a, b, c. To simplify imagine the equation is:

diff(y, t, 3) + a.*diff(y, t, 2) + b.*diff(y, t) + y.^(3/2) + c.*y == 0;

I know the discrete values of y and the time vector. I want to know tha values of a, b, c that solve this equation. Clearly there will be more combination of this values, but is there a way to solve it? (like with ode function or others).

One idea would be to express all the derivatives as function of y (like using finite differences) and then solve, is there a way to do it automatically on matlab?

[Sam Chak]

Hi @Raffaele

Can you confirm if the data pair is absolutely governed by the dynamics?

[Raffaele]

Actually the dynamic is more complicated, because there are also multiplications and divisions of the unknown parameters and also combination with other known parameters (like mass, spring stiffness, ecc.), like shown in my answer below. I just wrote that equation to simplify the presentation of the problem.

Answer 1

[Torsten]

Use "gradient" three times to get approximations for y', y'' and y''' in the points of your t-vector.

Let diff_y, diff_y2, diff_y3 be these approximations as column vectors of length n.

Now form a matrix A as

A = [diff_y2,diff_y,y]

and a column vector v as

v = [-y.^(3/2) - diff_y3]

Then approximations for a, b and c can be obtained using

sol = A\v

with

a = sol(1)
b = sol(2)
c = sol(3)

[Raffaele]

I'm sorry, maybe I should have been more precise in explaining the problem. Unfortunately I don't have the three unknown alone as in the equation above, but the coefficients are combination of unknown and known parameters.

diff(y, t, 3) + (m + 1/a).*diff(y, t, 2) + (k + (m*a) + R.*(b + c).*y.^(1/2)).*diff(y, t) + ...
    ((8*R^(1/2)*b)./(3*m*a)) .* y.^(3/2) + k.*y./(m*a) + (k/(m*a) + 1/a) == 0

This is how actually the equation looks like. Where a, b, c are the unknown and m, k and R are known.

[Sam Chak]

Hi @Raffaele

I understand that you want to fit the ODE into the data. Can you try out @Torsten's proposed method on this toy problem? See how correctly you can identify the 3 coefficients using sol = A\v.

% Coefficients
a = 60;
b = 1200;
c = 8000;

% System
f = @(t, y) [y(2); 
             y(3); 
           - a*y(3) - b*y(2) - c*y(1) - y(1)^(3/2)]; 

% Settings
y0 = [1 0 0];
tspan = [0 1];
[t, y] = ode45(f, tspan, y0);

plot(t, y(:,1)); grid on
ylabel('y(t)');
xlabel('Time, t');

[Torsten]

but the coefficients are combination of unknown and known parameters.

Then determine y',y'' and y''' as described and use "lsqnonlin" to determine the parameters.

As functions f_i for lsqnonlin, you have to specify the n function

f_i(a,b,c) = diff(y, t, 3) + (m + 1/a).*diff(y, t, 2) + (k + (m*a) + R.*(b + c).*y.^(1/2)).*diff(y, t) + ...
             ((8*R^(1/2)*b)./(3*m*a)) .* y.^(3/2) + k.*y./(m*a) + (k/(m*a) + 1/a)

[Raffaele]

@Torsten So you are suggesting to use lsqnonlin to find a, b, c, so give an initial guess for this three parameters, then find the solution in terms of y and minimize the error of this y simulated with the acquired one?

Because if I use initial guess for a, b, c, so that I know everything except y and its derivative, and then I use dsolve to find y(t) (to compare it to the acquired one), dsolve cannot find a solution to the problem.

[Raffaele]

@Sam Chak Mh I don't understand, since the coefficients are combination of known and unknown parameters, how can extract the unknown parameters using sol = A\v?

[Torsten]

Because if I use initial guess for a, b, c, so that I know everything except y and its derivative, and then I use dsolve to find y(t) (to compare it to the acquired one), dsolve cannot find a solution to the problem.

Your differential equation is nonlinear - thus "dsolve" will not succeed to solve it. Depending on the initial or boundary conditions, use "ode45" or "bvp4c".

[Raffaele]

And how can I solve an equation like that with ode45, considering there is a 3rd order derivative and also y^(3/2) and y^(1/2)?

[Torsten]

What are your 3 initial or boundary conditions ?

[Raffaele]

Since I know y(t), I know the initial values of y and its derivative. It should be y(0) = 0; y'(0) = -2.5; y''(0) = 100. But since there is y^(1/2) and y^(3/2), I thought to shift y, so using y(0) = 10. In any case I know the initial values but I have problems setting the ode solution.

[Torsten]

a = 2;
b = 3;
c = 4;
fun = @(t,y)[y(2);y(3);-a*y(3)-b*y(2)-y(1)^1.5-c*y(1)];
y0 = [0 ;-2.5 ;100];
tspan = [0 1];
[T,Y] = ode15s(fun,tspan,y0);
plot(T,Y)

Warning: Imaginary parts of complex X and/or Y arguments ignored.

grid on

[Raffaele]

Thank you, so now I try to insert it in an lsqnonlin to find the best option for a, b, c

[Torsten]

For the lsqnonlin fitting procedure, you can try two ways:

The first is without using the ode integrator and generating y',y'' and y''' from your data by using the "gradient" function.

The second is with the ode integrator by integrating the differential equation for given parameters a, b and c and by comparing the result with your given data for y.