Solving an initial value problem when one of the initial values is unknown

I am trying to numerically solve the ODE given by a*y"(t) + b*y'(t) = c, where a, b, and c are positive constants. I set this up as an initial value problem:

c(0) = c = 0.1

y(0) = 0

y'(0) = unknown, but I estimate it to be 0.1

Also:

dc/dt = 0

dy/dt = y(3)

d^2 y / dt^2 = (c-b*y(3))/a

This give me:

a = 10000;
b = 10000;
c = 0.1;
ode = @(t, y) [0; y(3); (c-b*y(3))/a];
y_dot_0_estimated = 0.1;
sol = ode45(ode, [0 1000], [c; 0; y_dot_0_estimated]);

My question is: is it possible to numerically figure out what the value for y'(0) should be? The solution seems to sensitive to this value, but I don't know what it is a priori. I know y(t) should approach c=0.1 at large t given large a and large b as given in the example, but this behavior might be different for other values of a or b.

Answer 1

[Torsten]

A second-order differential equation needs two boundary condition for y, either both at t = tstart (initial value problem) or as boundary conditions at t = tstart and t=tend (boundary value problem).

[L'O.G.]

@Torsten Both what? Do you mean both y(0) and y'(0)? Is there any way of numerically estimating y'(0) in my case?

[Torsten]

Initial value problem:

y(tstart) and y'(tstart) are given.

Boundary value problem:

y(tstart) and y(tend) are given or

y(tstart) and y'(tend) are given or ...

Answer 2

[John D'Errico]

Simple enough, since this ODE has an analytical solution, even without any boundary conditions specified. (If the ODE did not have an analytical soution, then it would still be possible to solve the problem, now using a tool like fzero, in combination with ODE45, or even lsqcurvefit. But the problem is simple enough to solve directly, so do as I show here.)

syms y(t)
a = 10000;
b = 10000;
c = 0.1;
dy = diff(y);
ddy = diff(dy);
y1 = dsolve(a*ddy + b*dy == c,y(0) == 0)

y1 = 

There is an ubdetermined constant, normally determined by the second boundary condition.

Now, suppose we assume the second BC is given as 0.1+delta, where delta is some perturbation. Now we can see how sensitive the result is to delta.

syms delta
y1prime_0 = subs(diff(y1,t),t,0)

y1prime_0 = 

So the first derivative of y1, at t==0 is just C_1, which is then just 0.1+delta. Is the result sensitvie to the value of delta? It does not seem like it should be that terribly sensitive at all.

syms C1
ysol = subs(y1,C1,0.1+delta)

ysol = 

That is the general solution of your problem, assuming you are uncertain about the exact first derivative BC at t==0.

So now lets look at the solution, when delta==0.

fplot(subs(ysol,delta,0),[0,1000],'b')

And I think now we see the problem. Do you see this problem as posed has a fast short term transient right near t==0?

hold on
fplot(subs(ysol,delta,0.01),[0,1000],'r')
fplot(subs(ysol,delta,-0.01),[0,1000],'g')
hold off

So a 10% perturbation to the value of y'(0) either way can cause a significant difference. The general curve has the same shape though. If we look at the short time behavior, the curves are:

fplot(subs(ysol,delta,0),[0,10],'b')
hold on
fplot(subs(ysol,delta,0.01),[0,10],'r')
fplot(subs(ysol,delta,-0.01),[0,10],'g')
hold off

Finally, can you "estimate" y'(0)? Yes, you could . As I showed, if you know the long term behavior of the curve, it would directly imply the value of y'(0). You can see the curve becomes linear in time for large t. So all you need to do is detemine the parameters of the straight line you should get for the long term behavior, that directly implies y'(0).

If all you had was some information on the uncertainty in y'(0). So, perehaps you decided that y'(0) has a mean of 0.1, but there was a standard deviation around that number, this too is doable. The result would be a distribution of solutions.

Anyway as I said, simple, almost trivial.