0

I'm studying numerical analysis and the answer for this problem (the question is in English) that is solved with the following code:

p1 = -20.2090;
p2 = 17.3368;
p3 = 272.9057;
beta = 7.8*pi/180;
x0 = 2000;

c = @(z)4800+p1+p2*z(1)/1000+p3*exp(-.75*z(1)/1000);
cp = @(z)(p2/1000-.75/1000*p3*exp(-.75*z(1)/1000));
q0 = (c(x0)/cos(beta))^2;

ode = @(t,z)[z(2); -q0*cp(z)/c(z).^3];
IC = [x0; tan(beta)];
[x,z] = ode45(ode,[0 25*6076], IC); 
plot(x,z(:,1))

But I don't understand the strategy of the above code. I could think that a second-degree differential equation should be solved using a transformation from second degree to a system of first-degree equation but it doesn't look like that is done. I do understand how we get the coefficients p1, p2, p3, beta and x0, that the c=... is the speed of sound related to depth and cp=... is the derivative. I would like to know what happens next and why the code solves the equation when it doesn't do a transformation. What is the purpose and meaning of z(2)in the ode and why does the ode look like is does? I know that ; in matlab means a second row and that z(2) means the second element in z, but why are those uder and how does it contribute to the correct answer?

Please help me understand.

Niklas Rosencrantz
  • 2,265
  • 1
    Frankly speaking, your doubts about transforming second order ODE to first order system weren't true. Just take a look at this line: ode = @(t,z)[z(2); -q0*cp(z)/c(z).^3] . It's exactly a transformation from ODE to system :) – Evgeny Sep 16 '15 at 12:23
  • 2
    The vector $z$ in your code has two components, which represent the value of the solution (call it $y$ for clarity) and the derivative of the solution, $y'$. Your ODE function says to use $y'$ (the second component) as the derivative of $y$, and to use the equation you were given for $y''$ to furnish the derivative of $y'$. You will always proceed this way when time-stepping an ODE of order higher than $1$, at least when using a "black box" function like ode45. – Ian Sep 16 '15 at 12:24
  • 1
    The link for "This problem" goes to a PDF. Do you want to tell us how to find that question, or were you expecting people to dig through that paper to find the question? – Thomas Andrews Sep 16 '15 at 12:24
  • 1
    @ThomasAndrews For me it opened at particular page with similar equations and functions. – Evgeny Sep 16 '15 at 12:25
  • 1
    The PDF has more than two pages. @Ian I count about eight. – Thomas Andrews Sep 16 '15 at 12:26
  • 1
    @ThomasAndrews You're right, I didn't catch that. But the link is supposed to redirect you to page 9, which is the only page which is mostly in English (except for its title). – Ian Sep 16 '15 at 12:27
  • 1
    Not every browser, apparently, supports such internal links to PDFs. @Ian – Thomas Andrews Sep 16 '15 at 12:28
  • Is it correct to say that the expression ode is the system matrix? Thank you for the comments and explanations. I think that I understand what the code does. It seems that ode45 is time-stepping with out system matrix and initial condition as arguments to ode45. Please correct me if I'm mistaken. – Niklas Rosencrantz Sep 16 '15 at 18:45
  • 1
    @Programmer400 It's not exactly a matrix; the function "ode" passed to ode45 is the function $f$ such that $y'=f(t,y)$ (where $y$ in this context is a vector). Your equation appears to be nonlinear, so $f$ is not given by a matrix. – Ian Sep 17 '15 at 00:13
  • Thank you for all the comments. You may take one of your comments and make it an answer. I feel that I understand much more now. – Niklas Rosencrantz Sep 17 '15 at 05:39

1 Answers1

1

Copying my previous comment:

The vector $z$ in your code has two components, which represent the value of the solution (call it $y$ for clarity) and the derivative of the solution, $y'$. Your ODE function says to use $y'$ (the second component) as the derivative of $y$, and to use the equation you were given for $y''$ to furnish the derivative of $y'$. You will always proceed this way when time-stepping an ODE of order higher than 1, at least when using a "black box" function like ode45.

Ian
  • 101,645