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

4 次查看（过去 30 天）
L'O.G.
L'O.G. 2023-10-11
编辑： Torsten 2023-10-12
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.

请先登录，再回答此问题。

回答（2 个）

Torsten
Torsten 2023-10-11
移动：Torsten 2023-10-11
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).
  2 个评论
L'O.G.
L'O.G. 2023-10-12
@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
Torsten 2023-10-12
编辑：Torsten 2023-10-12
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 ...

请先登录，再进行评论。

John D'Errico
John D'Errico 2023-10-12
编辑：John D'Errico 2023-10-12
Ran in:
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.

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

产品


版本

R2021b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by