How to solve an ODE with 2 initial conditions and 2 unknown parameters and 3 boundary conditions (overdetermined?)
8 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
e_frog
am 27 Okt. 2020
Kommentiert: e_frog
am 29 Okt. 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 Kommentare
Walter Roberson
am 28 Okt. 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
am 29 Okt. 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.
Akzeptierte Antwort
Alan Stevens
am 28 Okt. 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 Kommentare
Alan Stevens
am 28 Okt. 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.
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Ordinary Differential Equations finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!