You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
How to solve an ODE with 2 initial conditions and 2 unknown parameters and 3 boundary conditions (overdetermined?)
15 views (last 30 days)
Show older comments
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


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
:







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
Walter Roberson
on 27 Oct 2020
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
on 27 Oct 2020
I suspect "broken" means discontinuous, but it would help if we knew what m, F(t), tmax, K and v were!
Walter Roberson
on 27 Oct 2020
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
on 27 Oct 2020
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
on 27 Oct 2020
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?
Alan Stevens
on 28 Oct 2020
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?
Walter Roberson
on 28 Oct 2020
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)+12000000/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)+12000000/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)+12000000/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)))+198780000000 ...
/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
on 28 Oct 2020
@ 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
on 28 Oct 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
on 29 Oct 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.
Accepted Answer
Alan Stevens
on 28 Oct 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 Comments
e_frog
on 28 Oct 2020
ty, for the answer! could you please explain what the function FF = search does?
Alan Stevens
on 28 Oct 2020
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
on 28 Oct 2020
Why are there initial guesses for a and b in line 2 given? Where do those come from?
Alan Stevens
on 28 Oct 2020
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?
Alan Stevens
on 28 Oct 2020
Edited: Alan Stevens
on 28 Oct 2020
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
on 28 Oct 2020
Edited: e_frog
on 28 Oct 2020
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
on 28 Oct 2020
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
on 28 Oct 2020
what would be in those 2 last lines? The ODEs dont have a term with a third derivative term in them.
Alan Stevens
on 28 Oct 2020
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
on 28 Oct 2020
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
on 28 Oct 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.
e_frog
on 29 Oct 2020
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: 


More Answers (0)
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list:
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)