How to solve an ODE with 2 initial conditions and 2 unknown parameters and 3 boundary conditions (overdetermined?)

I have an ODE of the form

where a and b are unknown parameters, and F(t) is a fractional function: e.g.

The initial conditions are

and

Because a and b are unknown, the soultion of the ODE is depeding on a and b.

The solution x has to fulfill these boundary conditions in the time interval :

(where K is a known constant, e.g. )

(where )

( the maximum of the derivative of the solution is v, where v is a known velocity, e.g. )

The first problem is, that the linear system of equations to solve for a and b is overdetermined, because there are 3 eqantions for 2 unknowns. The second problem is, that i cant use the symbolic toolbox to compute the solution of the ODE, because of the broken rational function F(t) (it takes far too long, approx. more than 10 hours).

Would bvp4c or bvp5c be suitable for my problem? If so, then how to use them. And if not, is there anything else i could try?

[Walter Roberson]

What is a "broken rational function" ? It seems to have to do with the ratio of two polynomials, but I do not find information about the "broken" part ?

[Alan Stevens]

I suspect "broken" means discontinuous, but it would help if we knew what m, F(t), tmax, K and v were!

[e_frog]

sorry, that was lost in translation. I meant a fraction. I will edit the question!

[Walter Roberson]

If you can guarantee that the rational fraction F(t) is continuous over the t range of interest, then I see no reason not to use bvp4c or bvp5c

https://www.mathworks.com/help/matlab/ref/bvp5c.html#mw_a00814aa-1df0-4bd3-a378-d416355f6599 contains some hints that are intended to be guideance about which to use, but to be honest I do not know enough about either one of them for the hints to be meaningful to me at this time.

[e_frog]

thx for your help so far! If I am trying to solve the problem as a BVP I still dont know how to incorporate the condition , because in does not lay on the boundaries of my interval , but inbetween.

[Alan Stevens]

You have defined F in terms of x. x starts at 0, but you want it to end at 5. What happens at the point where x =1?

Or did you mean to define F(t) = (t+1)^2/(t-1)^2, in which case tmax is presumably less than 1?

[e_frog]

yes thank you, i meant t instead of x in F(t). I corrected the question

[Alan Stevens]

I've tried using fminsearch with various values of tmax < 1, but can't find values of a and b that achieve both x(t=tmax) = 5 and vmax = 3 at the same time.

You give K = 5 and vmax = 3 as examples, but are they the actual values you are interested in?

[e_frog]

Hi thanks for your persistance! The actual values would be:

[Walter Roberson]

If you use the initial equations

[ x(1/4) == 1/2, subs(diff(x(t), t), t, 1/4) == 0, 30000*diff(x(t), t, t) == a*diff(x(t), t) + b*x(t) + (111000*(t - 33/1000)^2)/(t + 37/1000)^2 + 882900]

then Maple says that the system has a closed form solution:

Ei = @expint
x(t) = 15664320000000*(b*exp(-1/60000000*(1000*t+37)*(-a+(a^2+120000*b)^(1/2))) ...
*(a-(a^2+120000*b)^(1/2)[phone]/7)*Ei(1,-1/60000000*(1000*t+37)*(-a+(a^2+     ...
120000*b)^(1/2)))-Ei(1,287/60000000*a-287/60000000*(a^2+120000*b)^(1/2))*(a-(a^ ...
2+120000*b)^(1/2)[phone]/7)*b*exp(-1/60000000*(1000*t+37)*(-a+(a^2+120000*b)^ ...
(1/2)))+((-50000/1813*b-99390000000/1813)*(a^2+120000*b)^(1/2)+(-50000/1813*a+  ...
60000000/287)*b-99390000000/1813*a)*exp(-1/240000*(4*t-1)*(-a+(a^2+120000*b)^(1 ...
/2)))-(a+(a^2+120000*b)^(1/2)[phone]/7)*exp(1/60000000*(1000*t+37)*(a+(a^2+   ...
120000*b)^(1/2)))*b*Ei(1,1/60000000*(1000*t+37)*(a+(a^2+120000*b)^(1/2)))+Ei(1, ...
287/60000000*a+287/60000000*(a^2+120000*b)^(1/2))*b*(a+(a^2+120000*b)^(1/2)+    ...
12000000/7)*exp(1/60000000*(1000*t+37)*(a+(a^2+120000*b)^(1/2)))+((-50000/1813* ...
b-99390000000/1813)*(a^2+120000*b)^(1/2)+(50000/1813*a-60000000/287)*b+         ...
99390000000/1813*a)*exp(1/240000*(4*t-1)*(a+(a^2+120000*b)^(1/2)))[phone] ...
/1813*(a^2+120000*b)^(1/2))/(a^2+120000*b)^(1/2)*b^2/(a-(a^2+120000*b)^(1/2))^3 ...
/(a+(a^2+120000*b)^(1/2))^3

This does not incorporate the condition

You can differentiate twice and try to solve for 0, but that is going to be tough going. And you have three parameters to the equation (a, b, t) so if you are talking about the maximum of the derivative, you have too few equations to solve exactly, so you are going to have to do a maximization that is constrained in t but unconstrained in a and b.

My numeric tests suggest strongly that the maximum of the derivative occurs at t = tmax . However, the best value I have been able to find so far is a maximum on the order of 1e-11, not even close to 3.6, so I am not at all certain the 3.6 can be achieved.

Do you have ideas on the maximum and minimum reasonable values for a and b ?

[Alan Stevens]

@ Walter The values my program produced are a = 1.521*10^5 and b = -5.8*10^-6. These produced the graph shown below.

[Walter Roberson]

@Alan: if Maple is correct about the explicit solution to the ode, then the locations you show are not even close to the values you show in the graph... Not unless we are working with completely different constants. The poster changed some of the values a lot along the way and I was using the constants from https://www.mathworks.com/matlabcentral/answers/627358-how-to-solve-an-ode-with-2-initial-conditions-and-2-unknown-parameters-and-3-boundary-conditions-ov#comment_1092173

Mind you, my tests did not have access to the information that a and b are positive; the actual fitting may be significantly worse than my earlier efforts.

[Alan Stevens]

@Walter: My values for a and b are quite possibly wrong - they only get the maximum xdot reasonably (Looks like I was trying for xf = 0.5 instead of 5. Hmm, I can't even see where the value of 0.5 was specified now!). As you say, the spec changed a few times. No matter what initial guesses I made for a and b, fminsearch always returned positive a and negative b, even if my initial guesses were both positive.

Answer 1

[Alan Stevens]

The following is as close as I can get. You might be able to do better by altering tolerances using setoptions, but I'll leave that to you:

tspan = [0 0.25];
a = 1.5*10^5; b = -5*10^6;
AB0 = [a; b];
AB = fminsearch(@search,AB0,[],tspan);

a = AB(1); b = AB(2);
disp('    a          b')
disp([a, b])

IC = [0 0];
[t, X] = ode45(@(t,x) fn(t,x,AB),tspan,IC);
x = X(:,1);
v = X(:,2);
 subplot(2,1,1)
 plot(t,x),grid
 xlabel('t'),ylabel('x')
 subplot(2,1,2)
 plot(t,v),grid
 xlabel('t'),ylabel('xdot')


function FF = search(AB,tspan)
         K = 0.5;
         vmax = 3.6;
         [~,X] = odesol(AB,tspan);
         xf = X(end,1);
         vm = max(X(:,2));
         vf = X(end,2);
         
         FF = (K-xf)^2 + (vmax - vm)^2 + vf^2;
end


function [t,X] = odesol(AB,tspan)
        IC = [0 0];
        [t, X] = ode45(@(t,x) fn(t,x,AB),tspan,IC);
end

function dXdt = fn(t,X,AB)
         m = 30000;
         a = AB(1);
         b = AB(2);
         x = X(1);
         v = X(2);
         dXdt = [v;
                 (a*v + b*x + 1110000*(t-0.0331)^2/(t+0.0371)^2 + 882900)/m];
end

This produces

[e_frog]

ty, for the answer! could you please explain what the function FF = search does?

[Alan Stevens]

This calculates the values of the parameters of interest, namely xf = x(t = 0.25); vmax = maximum(xdot) and vf = v(t = 0.25), for given values of 'a' and 'b' as supplied by fminsearch.

It does so by calling odesol, which, in turn uses ode45 (in fact, I could have missed out odesol and put the ode45 function directly into the search function!). The function search then forms the difference between the desired values of xf, vmax and vf and their computed values, squares the differnces and adds them up, and returns the result to fminsearch.

Fminsearch then adjusts 'a' and 'b' in an effort to make the returned sum of squares as close to zero as it can (If the returned sum of squares were to be exactly zero, the calculated parameters would exactly match the desired values.

[e_frog]

Why are there initial guesses for a and b in line 2 given? Where do those come from?

[Alan Stevens]

fminsearch requires initial guesses for the parameters it's trying to optimise. I just did a little numerical experimentation (i.e. guesswork!) before settling on those. Feel free to try other values.

Do you have an idea of what the final values of a an b should be?

[e_frog]

a and b must be positive. That's the only condition for both values

[Alan Stevens]

I've tried with both a and b having positive initial guesses, and the routine always comes up with the final values of 'a' positive and 'b' negative.

Double check the equations. The forcing function F(t) is always positive and has some large values, so one would expect at least one of a and b to be negative in order to stop a runaway acceleration.

[e_frog]

Thx, thats what I thought too, solving the ODE numerically! I still need some help with the function

dXdt = fn(t,X,AB).

How would I modify it for calculating a system of ODEs?

I tried the following;

function dXdt = fn(t,X,AB)
         [A1,B1,C1,D1] = somefunction();
         m1 = 30000;
         m2 = 300;
         d = AB(1);
         c1 = AB(2);
         c2 = AB(3);
         s1 = X(1);
         s1_dot = X(2);
         s2 = X(3);
         s2_dot = X(4);
         dXdt = [s1_dot;
                 -(d/m1)*s1_dot-(c1/m1)*s1*(d/m1)*s2_dot+(c1/m1)*s2+(1/m1)*(-A1*((t-B1)^2)/((t+C1)^2)+D1);
                 s2_dot;
                 (d/m2)*s1_dot+(c1/m2)*s1-(d/m2)*s2_dot-((c1+c2)/m2)*s2];
end

[Alan Stevens]

function dXdt = fn(t,X,AB)
         [A1,B1,C1,D1] = somefunction();  % What is this? A1 etc don't seem to be used below!
         m1 = 30000;
         m2 = 300;
         d = AB(1);
         c1 = AB(2);
         c2 = AB(3);
         s1 = X(1);
         s1_dot = X(2);
         s2 = X(3);
         s2_dot = X(4);
         dXdt = [s1_dot;
                 -(d/m1)*s1_dot-(c1/m1)*s1*(d/m1)*s2_dot+(c1/m1)*s2+(1/m1)*(-A1*((t-B1)^2)/((t+C1)^2)+D1);
                 s2_dot;
                 (d/m2)*s1_dot+(c1/m2)*s1-(d/m2)*s2_dot-((c1+c2)/m2)*s2;
                  this line will contain the expression for ds1_dot/dt;
                  this line will contain the expression for ds2_dot/dt];
end

[e_frog]

what would be in those 2 last lines? The ODEs dont have a term with a third derivative term in them.

[Alan Stevens]

You have extracted four values from X. That means four gradients must be returned. I assumed you had equations for the rate of changes with time of s1_dot and s2-dot. In general if you have a system of first order equations (which you must have in order to use ode45 - they cannot be higher order) you must return the same number of gradients as the number of terms in X.

[e_frog]

I actually am not sure, what i am doing there with X. I just tried to replicate what you did for one equation. The two ODE are both second order ODEs

[Alan Stevens]

Originally, you had one second order ode which was turned into two first order odes - necessary before being offered to ode45.

If you now have two second order odes you must turn them into four first order odes, then you can offer them to ode45. If your two 2nd order odes are, say d2x/dt2 = f(...) and d2y/dt2 = g(... ), then you will need to make them

dx/dt = v; dv/dt = f(...); dy/dt = w; dw/dt = g(...)

The values that get passed to the gradient function (which in my case I just called fn) will be X = [x v y w]

The values that are returned will be dXdt = [v; f(...), w, g(...)] - (note, this must be a column vector).

Hope this helps.

[e_frog]

I made both the version with 1 ODE and 2 ODEs work for me. But the downside is, that the solutions never match the boundary conditions exactly.

For example: for 1 ODE the boundary condition is not met, the resut for that is always how do I adjust the tolerances so that the conditions are met better (or maybe met exactly)? For the system with 2 ODEs the conditions are way off.

E.g.: Goal: , but actual result: