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일
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 CenterFile Exchange에서 Numeric Solvers에 대해 자세히 알아보기

태그

제품


릴리스

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by