Using ode45 for a non linear differential equation with a term in a vector form

Hello,

I want to solve the 2nd order differential equation:

d^2(xus)/dt^2 = -K/Mus*xus-Ms/Mus*d^2(xs)/dt^2

I found as = d^2(xs)/dt^2 experimentally with potentiometer so I have as a 90000x1 vector that corresponds to the acceleration at each time t and now I want to find xus in function of this same time t using ode45

The problem is that I don't know how to deal with as in ode45 and I have the error:

Unable to perform assignment because the size of the left side is 1-by-1 and the size

of the right side is 0-by-1.

Considering that as is given the rest of my code is:

Mus=100;

Ms=10;

K=1000;

tspan=linspace(0,(size(as)-1)*0.02,size(as));

x0 = [0 0];

[t,x] = ode45(@(t,x) fun(t,x,as,Mus,Ms,K,tspan), tspan, x0);

function dx_dt = fun(t,x,as,Mus,Ms,K,tspan)

i=find(tspan==t);

a=as(i);

dx_dt(1,1) = x(2);

dx_dt(2,1) =-K/Mus*x(1)-Ms/Mus*a;

end

Thank you in advance for your help

Answer 1

[Star Strider]

I do not understand the problem.

Is the objective to estimate a parameters of the differential equation from data?

The time value at each iteration is supplied to the ode45 function so there is no reason to specify it directly. If the data were collected at specific times, use that time vector for ‘tspan’ in the differential equation call. Are the values of ‘K’, ‘Ms’, and ‘Mus’ known, or are they to be estimated (as parameters) from the data?

The derivation is fairly straightforward using the Symbolic Math Toolbox —

syms D2xs K Ms Mus T t xus(t) Y 
Dxus = diff(xus);
D2xus = diff(Dxus);
Eqn = D2xus == -K/Mus*xus-Ms/Mus*D2xs

Eqn(t) = 

[VF,Sbs] = odeToVectorField(Eqn)

VF = 

Sbs = 

xusfcn = matlabFunction(VF, 'Vars',{T,Y,D2xs,K,Ms,Mus})

xusfcn = function_handle with value:
    @(T,Y,D2xs,K,Ms,Mus)[Y(2);-(K.*Y(1))./Mus-(D2xs.*Ms)./Mus]

It would help to upload/attach the time and vectors (as a matrix as (90000x2) matrix if the intent is to estimate parameters from it.

.

[Maxime Jean Marc Lenoir]

Thank you for your answer.

The objective is to find 'xus' according to the differential equation.

‘K’, ‘Ms’, and ‘Mus are known parameters and 'as' is determined experimentally with a samplig period T=0.02 seconds.

Therefore, I would like to obtain a vector 9000x1 for 'xus' that satisfies the differential equation at each sampling time.

[Star Strider]

I am a bit lost.

If the data are , all the parameter values are known. and the objective is to find , I would inmtegrate it 2 times with trapz and use that. Plot that result (perhaps estimating constants-of-integration as well at eash step, or estimating them using appropriate linear or nonlinear regression techniques) with the computed result of the numerically-integrated differential equation.

If the objective is to estimate the unknown values of ‘K’, ‘Ms’, and ‘Mus’, that is different and requires nonlinear parameter estimation of the differential equation as it is being integrated. Both of these problems are relatively straightforward to solve, however require significantly different approaches.

What is the desired outcome?

.

[Maxime Jean Marc Lenoir]

Let me redefine the problem, the differential equation is:

K, Ms and Mus are given parameters in the problem, is found experimentally with a sampling period of T=0.02s and its values at each time t=nT with n=0,1,...,8999 are stored in the 9000x1 vector 'as'.

Finally I want to find in a vector form with the values at each time t=nT.

I thought of using ode45 to solve the differential equation and find but there are some mistakes in the code since I don't know how to deal with the vector 'as' in ode45.

[Star Strider]

So there is nothing to be estimated, all the parameters are known, and the desired result is to calculate as a function of vector if I understand correctly.

The differential equation can be solved analytically —

syms D2xs K Ms Mus T t xus(t) Y xus0 Dxus0 
sympref('AbbreviateOutput',false);
Dxus = diff(xus);
D2xus = diff(Dxus);
Eqn = D2xus == -K/Mus*xus-Ms/Mus*D2xs

Eqn(t) = 

cndx = [xus(0)==xus0, Dxus(0)==Dxus0];
xus = dsolve(Eqn, cndx)

xus = 

xus = simplify(xus, 500)

xus = 

xusfcn = matlabFunction(xus)

xusfcn = function_handle with value:
    @(D2xs,Dxus0,K,Ms,Mus,t,xus0)(exp((t.*sqrt(-K.*Mus))./Mus).*(K.*xus0+D2xs.*Ms+Dxus0.*K.*Mus.*1.0./sqrt(-K.*Mus)))./(K.*2.0)+(exp(-(t.*sqrt(-K.*Mus))./Mus).*(K.*xus0+D2xs.*Ms-Dxus0.*K.*Mus.*1.0./sqrt(-K.*Mus)))./(K.*2.0)-(D2xs.*Ms)./K

From here, supply the necessary parameter values, the ‘t’ and (denoted here as ‘D2xs’) values to get the result.

Note that the ‘xus0’ and ‘Dxus0’ parameters are the initial conditions for the differential equation, so set them to whatever they are supposed to be. They must be supplied. Neither appear in the denominator of any fraction, so it is highly unlikely that there will be ±Inf or NaN results if they are set to 0.

In this instance, it is not necessary to do a numerical integration.

.

[Maxime Jean Marc Lenoir]

Great ! Thank you so much for your help

[Star Strider]

As always, my pleasure!

The order of the arguments can be changed by using the 'Vars' name-value pair in the matlabFunction call. That might make them easier to work with. (All calls to the function have to have their argument list changed to match it, if it the order is changed.)

.