MATLAB Answers

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

105 views (last 30 days)
e_frog
e_frog on 27 Oct 2020 at 7:30
Commented: e_frog on 29 Oct 2020 at 14:08
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?

  13 Comments

Alan Stevens
Alan Stevens on 28 Oct 2020 at 10:47
@ 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
Walter Roberson on 28 Oct 2020 at 23:39
@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
Alan Stevens on 29 Oct 2020 at 8:33
@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.

Sign in to comment.

Accepted Answer

Alan Stevens
Alan Stevens on 28 Oct 2020 at 9:36
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

  13 Comments

e_frog
e_frog on 28 Oct 2020 at 18:40
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
Alan Stevens on 28 Oct 2020 at 19:21
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
e_frog on 29 Oct 2020 at 14:08
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:

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by