Initial conditions for solving numerically a second order ODE

Hello,

I am trying to solve numerically a second-order ODE and I want to implement the following intial conditions:

y(0)=0 and y'(0)=-1. How can I implement those initial conditions? The code and the ODE are given below. MAong others I already gave a try by importing initial condition cond=Dy(0)==-1 and yInit=[0 0] but it seems it did not work.

clear all
syms y(t)
Dy = diff(y);
eqn = diff(y,2)==-104.17444*diff(y)-94.43862*((1/((1+y/0.000878)^0.5))-(1/((1+y/0.000878)^2)))-9.81+0.0002*exp(500*t);
cond=Dy(0)==-1;
V = odeToVectorField(eqn,cond);
M = matlabFunction(V,'vars',{'t','Y'});
interval = [0 0.01];
yInit=[0 0];
ySol = ode45(M,interval,yInit);

Any suggestions would be very much appreciate it.

Thanks!

Best,

Christos

[darova]

How do you know it doesn't work? Did you try to plot the solution?

[Christos Stamatopoulos]

Dear darova,

Just to clarify, of course I get a numerical solution of the ODE. What is not working is that the initial conidtions that I want are not imlemented at the solution properly.

Best,

Christos

Answer 1

[John D'Errico]

Several problems in your code.

First, you don't pass the initial conditions into odeToVectorField. It is not solving your ODE for you, just converting it into a pair of differential equations, something that could have been done by hand. Regardless, you don't pass in the initial conditions.

Next, using tools like matlabFunction and odeToVectorfield does not teach you anything about what you doing, nor even make you think about it. Instead, you allow the computer to do the thinking for you. So learn how to solve the problem first, rather than relying on the computer. Otherwise, you will make lots of mistakes, because you are passing virtually random things into the tools you don't understand how to use.

You convert a second oder ODE into a first order system by a simple transformation. Thus, if you have the problem:

y'' = f(t,y,y')

then you create a second unknown variable, to replace y'. I suppose I could now have a problem with two unknowns, y and yp. Here yp is the new variable I have just created. both y and yp are a function of time.

y' = yp

yp' = f(t,y,yp)

You agree that what I have created is now a simple system of first order ODES? Looking at your problem:

eqn = diff(y,2)==-104.17444*diff(y)-94.43862*((1/((1+y/0.000878)^0.5))-(1/((1+y/0.000878)^2)))-9.81+0.0002*exp(500*t);

pretty(vpa(eqn,5))
  2
 d                                       94.439                   94.439                   d
--- y(t) == exp(500.0 t) 0.0002 + -------------------- - ----------------------- - 104.17 -- y(t) - 9.81
  2                                                  2   sqrt(1139.0 y(t) + 1.0)          dt
dt                                (1139.0 y(t) + 1.0)

First, I will argue that you very, very likely will have numerical problems if you go out too far in time. That exp(500*t) in there will explode, doing nasty things to any code. I suppose you may survive, as the maximum time is just 0.01, so we will have numbers as large as exp(5). Now, re-write your problem in terms of the new variable, as a pair of equations.

y' = yp

yp' = -104.17444*yp-94.43862*((1/((1+y/0.000878)^0.5))-(1/((1+y/0.000878)^2)))-9.81+0.0002*exp(500*t)

What does a numerical solver expect? It needs a function of the form f(t,Y), where Y is a vector of length 2 in this case. Think of y(1) as your original variable y, and Y(2) as the new variable yp. I'll create a function handle as:

FUN = @(t,Y) [Y(2); ...
       -104.17444*Y(2)-94.43862*((1/((1+Y(1)/0.000878)^0.5))-(1/((1+Y(1)/0.000878)^2)))-9.81+0.0002*exp(500*t)];

If you look carefully at what I did, I replaced diff(y) in your equation with the variable Y(2). y in your equation becomes Y(1).

Next, what are your initial conditions? They are not differential equations, even though you seem to be thinking of them as differential equations or at least treating them as such. They are initial conditions - thus a fixed vaue of your function (or its derivative) at a fixed point in time.

We had y(0) = 1, and Dy(0) = yp(0) = -1.

yinit = [0,-1];
tspan = [0,0.01];

I have no idea why you decided to set them both to 0 in your code, when you explicitly stated that dy(0) was -1.

[t,ySol] = ode45(FUN,tspan,yinit);

Looking at the result even for t only 0.01, we see:

ySol(end,:)
ans =
   3.4325e+12   1.7163e+15

That exponential probably started to drive things rapidly to hell. I was wondering if the problem would be perceived as a stiff one, but using a stiff solver did not change things a lot, nor did ODE45 hang up.

Is the result a meaningful one? Hard to know, but there are numerical risks in this. ySol seems to be diverging pretty fast here. I would certainly suspect the physical meaningfulness of a solution that starts out as smoothly as it does, but then diverges rapidly to infinity.

Finally, could you have done this using odeToVectorField and matlabFunction? Well, yes. But learn to solve the problem first before you just throw a problem you don't understand into some tool you don't understand either.