How to solve an ODE with 2 initial conditions and 2 unknown parameters and 3 boundary conditions (overdetermined?)
8 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
e_frog
el 27 de Oct. de 2020
Comentada: e_frog
el 29 de Oct. de 2020
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 comentarios
Walter Roberson
el 28 de Oct. de 2020
@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
el 29 de Oct. de 2020
@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.
Respuesta aceptada
Alan Stevens
el 28 de Oct. de 2020
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 comentarios
Alan Stevens
el 28 de Oct. de 2020
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.
Más respuestas (0)
Ver también
Categorías
Más información sobre Ordinary Differential Equations en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!