Solving first order differential equation

Hello, I've tried multiple times to solve the following differential equation in Matlab but no luck so far. I have about 131 different values of U for 131 seconds of time t. A, B, r are constants, y and dy/dt has initial conditions of 0. I want to calculate L for each time t and plot a graph.

I've tried using ode23 and created a function. But it doesn't work. Any help in coding this in Matlab would greatly appreciated as I'm getting desperate! Thank you very much.

[Stephan]

Please attach what you have tried

[Jon]

Before giving you advice advice about how to solve the problem using MATLAB I need to understand what is known and what is being solved for. You say the the load L is an output, so I assume that is unknown. Do you know U as a function of time? If so how is it specified. Do you have it defined by a function, or are you given values of U for specific values of time? You say that h is a displacement, so is the velocity U the rate of change (derivative with respect to time) of h?

Answer 1

[Star Strider]

You need to re-formulate your ‘first_order’ function. MATLAB derives a different set of equations:

syms A B h(t) L r U h0 Y 
Eqn = diff(h,2) == 2/(U*r*B)*L + U*A + 0.6211*U*h/B;
[VF,Sbs] = odeToVectorField(Eqn)
first_order = matlabFunction(VF, 'Vars',{t,Y,A,B,L,r,U})

producing:

VF =
                                          Y[2]
 A*U + (6211*U*Y[1])/(10000*B) + (2*L)/(B*U*r)

Sbs =
  h
 Dh

first_order = @(t,Y,A,B,L,r,U)[Y(2);A.*U+(U.*Y(1).*6.211e-1)./B+(L.*2.0)./(B.*U.*r)];

[Saruultugs Batzorig]

Thank you for your answer, do I also need need to change the simulate.m file to fit all the parameters as ode23 = (@(t,Y,A,B,L,r,U)first_order(t,Y,A,B,L,r,U), t, y0)? Also, I made a mistake in my original question, I put a second order I've changed it to first order.

[Star Strider]

My pleasure.

I do not understand what you intend by ‘fit all the parameters’. Do you have data that you want to fit your differential equation to in order to estimate parameters, or do you simply want to loop through the various values of the parameters that you have already chosen values for?

[Saruultugs Batzorig]

I have about 131 different values of U for 131 seconds of time t. A, B, r are constants, y and dy/dt has initial conditions of 0. I want to calculate L for each time t. I hope that clears things up.

[Star Strider]

I let the Symbolic Math Toolbox do the ‘heavy lifting’, since your differential equation turned out to have an analytic solution:

syms A B h(t) L r U h0 Y 
Eqn = diff(h) == 2/(U*r*B)*L + U*A + 0.6211*U*h/B;
hsol = dsolve(Eqn, h(0)==h0)                                        % Analytic Solution
Lsol = solve(hsol,L);                                               % Solve For ‘L’
Lsol = simplify(Lsol, 'Steps',500)                                  % Simplify
Lfcn = matlabFunction(Lsol)                                         % Function Handle

producing:

hsol =
-(20000*L + 10000*A*B*U^2*r - 6211*U^2*r*exp((6211*U*t)/(10000*B))*(h0 + (10000*A*B*r*U^2 + 20000*L)/(6211*U^2*r)))/(6211*U^2*r)

Lsol =
- (U^2*r*(6211*h0 + 10000*A*B))/20000 - (6211*U^2*h0*r)/(20000*(exp((6211*U*t)/(10000*B)) - 1))

Lfcn = @(A,B,U,h0,r,t)U.^2.*r.*(h0.*6.211e+3+A.*B.*1.0e+4).*(-5.0e-5)-(U.^2.*h0.*r.*3.1055e-1)./(exp((U.*t.*6.211e-1)./B)-1.0);

The ‘Lfcn’ anonymous function should give you ‘L’ as a funciton of the other variables.

I have no idea what your data are or what you are doing. If you have data for ‘h(t)’ as well, the derivation changes to:

Lsol = solve(h == hsol,L);                                          % Solve For ‘L’
Lsol = simplify(Lsol, 'Steps',500)                                  % Simplify
Lfcn = matlabFunction(Lsol)                                         % Function Handle

with appropriate changes in the functions to give:

Lfcn = @(A,B,U,h0,r,t)U.^2.*r.*(h0.*6.211e+3+A.*B.*1.0e+4).*(-5.0e-5)-(U.^2.*r.*(h0-h(t)).*3.1055e-1)./(exp((U.*t.*6.211e-1)./B)-1.0);

I will let you sort that.

[Saruultugs Batzorig]

That was really helpful, thank you so much.

[Star Strider]

As always, my pleasure!

[Saruultugs Batzorig]

One last question, if A were an integral given as:

how would I implement in the differential equation eqn? Here is my attempt but it doesn't seem to work. I used the trapz method with linspace values.

Eqn = diff(h) == 2/(*U*r*B)*L +...
    (trapz(exp((-U*t/B)*0.1843*0.3)*h,linspace(0,t)))/(-3.7549)) ...
    + 0.6211*U*h/B;

[Star Strider]

This has now become an integro-differential equation.

I have not worked with integro-differential equations in decades. However, it seems reasonable to code it as:

A = @(t) integral(@(t) exp((-U*t/B)*0.1843*0.3)*h, 0, t);

either inside your ‘first_order’ function, or calling it from within your ‘first_order’ funciton. In that instance, you may need to include the other variables as well:

A = @(t,U,B,h) integral(@(t) exp((-U*t/B)*0.1843*0.3)*h, 0, t);

Then call it appropriately.

Note: I did not test this, so I am posting it as UNTESTED CODE.

[Saruultugs Batzorig]

Do you mean inside the function as in like this? Sorry for all the questions, my Matlab knowledge is very basic.

syms A B h(t) L r U h0 Y
A = @(t,U,B,h) integral(@(t) exp((-U*t/B)*0.1843*0.3)*h, 0, t);
Eqn = diff(h) == 2/((-3.7549)*U*r*B)*L +U*((0.1843*U)/(B^2))*A/(-3.7549) ...
    + (0.6211/(-3.7549))*U*h/B;
hsol = dsolve(Eqn, h(0)==h0);                                  
Lsol = solve(hsol,L);                                         
Lsol = simplify(Lsol, 'Steps',500);                            
Lfcn = matlabFunction(Lsol);     

[Star Strider]

I was thinking of an end result that (with a bit of editing) would be:

A = @(t,U,B,h) integral(@(t) exp((-U*t/B)*0.1843*0.3)*h, 0, t);
Lfcn = @(B,U,h0,r,t)(A(t,U,B,h).*U.^3.*r.*(-9.215e-2))./B+(U.^2.*h0.*r.*3.1055e-1)./(exp((U.*t.*1.654105302404858e-1)./B)-1.0);

Note that ‘h’ is the left-hand-side of your original ‘hsol’ equation, so it should also appear in ‘Lsol’, although now in the RHS expression, and in the argument list of ‘Lfcn’.

[Saruultugs Batzorig]

Thank you! Is there a way to use Lfcn without having to input the value for h as I do not know, and must be found from the equation itself?

[Star Strider]

Is ‘L’ a function or a parameter? (I always considered it a parameter.) There may be ways to solve for ‘L’ as a parameter, even varying with respect to time, if ‘h’ and ‘t’ are otherwise specified. This then becomes a parameter estimation problem.

[Saruultugs Batzorig]

L is a parameter to be found from the equation yes.t is specified, unfortunately only the initial values for h is known and is 0.

[Star Strider]

Parameters by definition do not vary. If you specify a value for ‘h’ at a particular time ‘t’, you can likely solve for ‘L’ at those values of ‘h’ and ‘t’. Do you have values for ‘h’ and ‘t’ otherwise? (Note that in that event, ‘L’ becomes a variable, and not a parameter.)

[Saruultugs Batzorig]

I see, L is a variable in that case. t being time goes from 0 to 100. I do not have values for h.

[Star Strider]

To estimate ‘L’ as a function of ‘t’, you would need to provide values for ‘h(t)’. You could then solve for ‘L’.

Alternatively, you could re-formulate your differential equation in terms of . However this still leaves the problem of defining ‘h’, unless you reformulate it to eliminate ‘h’.

I have no idea what you are doing.

Answer 2

[Jon]

Your verbal problem description in which you call L and ouput, but in fact you know what it is, somewhat confuses me, but in anycase going from your equations I think I see what you are trying to do.

Besides the details of how to get extract the variables from what is returned by ode23, I think you have some more fundamental issues.

You are solving for dh/dt , but you want h. If you want h then you should set up a system of two first order ode's

So for example defining y1 = h, y2 = dh/dt

dydt = fun(t,y)
    dydt = [y(2)
        (2/(U.*r*B)).*L+U.*A+(U./B)*0.6211*y(1)];

Note that you can also use two left hand argument to ode23 and get the answers returned as y and t rather than in a structure if you prefer.

The additional arguments can be sent using the anonymous function as you did before as long as they are in the local workspace when you call ode23.

Sorry I see now that @Starstrider has been giving you additional ideas while I was writing this answer, but maybe this is helpful too.

[Saruultugs Batzorig]

Thank you for you answer, ultimately I'm trying to plot L against time t for the given equation with the said initial conditions.

[Jon]

What exactly do you know, and what are you trying to solve for?

You can not solve for both L and h knowing only U as a function of time.

If you know h as a function of time and U as a function of time you could differenitate h twice then substitute for d2hdt in your function and solve for L as a function of time. (Not solving an ODE at all).

If you know L(t) as a function of time and U(t) as a function of time then you could use ode23 but you must include them within your dydt function. So just as an illustrative made up example if U(t) = sin(t) and L(t) = cos(t) you would have

dydt = fun(t,y)
U = sin(t);
L = cos(t)
dydt = fun(t,y)
    dydt = [y(2)
        (2/(U.*r*B)).*L+U.*A+(U./B)*0.6211*y(1)];