**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?)

24 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

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

Walter Roberson
on 27 Oct 2020

Alan Stevens
on 27 Oct 2020

e_frog
on 27 Oct 2020

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

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

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 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

### 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

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

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

e_frog
on 28 Oct 2020

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

### 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.