How to Solve a Shooting Method Problem with Specific Boundary Conditions?

[electronic engineer]

Hello, can anyone give me the general instructions of solving shooting method problem:

dy1/dx=-y1^2*y2
dy2/dx=y1*y2^2

with the boundary conditions: y1(0)=1, y2(1)=2

 

[phyzguy]

To solve a set of differential equations like this numerically, you generally start at some point and then step forward through your domain of interest, using a method like Euler's method or Runge-Kutta. The problem here is that you don't know the starting boundary conditions at any point, you know them at two different points. So to use the shooting method, you start at some point, like x=0. You know y1(0) = 1, but you don't know y2(0). So you guess a value, say y2(0)=1. Then you numerically step forward to x=1, and see what the value of y2 is. If y2(1) =2, you are done, but probably it is greater or less than 2. Based on y2(1), you adjust your initial guess of y2(0) and step forward again. Hopefully this time you are closer. You keep iterating until y2(1) = 2 to the desired degree of precision. Alternatively, you could start at x=1, guess a value for y1(1) and step down to x=0. The process is the same.

 

[electronic engineer]

Thanks but still not clear how to step forward to the other boundary point. We have two equation, i don't know from which one to start and how. What assumptions should be made and what concrete processes? Can you make it a little bit more clear?! Thank you.

 

[phyzguy]

Are you familiar with Euler's method? It is well explained there, but I'll summarize. You know y1(0) and y2(0). From your Diff. Eq, you can calculate y1'(0) and y2'(0). Since y1' = Δy/Δx, you can write Δy1 = y1'(0) * Δx, Δy2 = y2'(0) * Δx. Then y1(Δx) = y1(0) + Δy1, y2(Δx) = y2(0) + Δy2. You just keep doing this, incrementing x by small steps (which I called Δx, and the Wiki article calls h), until you arrive at x=1.

Edit: Corrected the expressions for Δy1 and Δy2.

 

[Chestermiller]

Are you allowed to solve this problem analytically?

Chet

 

[phyzguy]

Are you allowed to solve this problem analytically?

Chet

A good question. I kind of got interested in this problem. It can be solved both numerically and analytically, and it is very instructive to compare the two methods and see how close the numerical method comes to the analytical solution.

 

[electronic engineer]

Are you allowed to solve this problem analytically?

Chet

No, only numerically.

 

[Chestermiller]

In iterating to get the initial value of y2 at x = 0, you are essentially solving a non-linear algebraic equation. What I usually do is apply Newton's method, approximating the derivative of the function numerically. Another approach would be to make use of the half-interval method (if you can get the solution for y2(0) bracketed).

Chet

 

[pyroknife]

You can find all sorts of shooting method BVPs online. It essentially turns a BVP into an initial condition problem

The problem you posted is a very simple one and tons of examples can be found on google.

 

[pasmith]

Thanks but still not clear how to step forward to the other boundary point. We have two equation, i don't know from which one to start and how. What assumptions should be made and what concrete processes? Can you make it a little bit more clear?! Thank you.

You have to treat them as the single vector equation $$\begin{pmatrix} y_1' \\ y_2' \end{pmatrix} = \begin{pmatrix} -y_1^2 y_2 \\ y_1y_2^2 \end{pmatrix}.$$ Every numerical ODE solution method deals with vector ODEs.

The point of the exercise is to take a guess for $y_2(0)$ and see whether $y_2(1) = 2$. You do need to do some analysis of the phase space first in order to obtain a plausible initial guess. Here the phase space is the $(y_1,y_2)$-plane. You can see that both the $y_1$-axis and the $y_2$-axis consist of fixed points, so a trajectory which starts with $y_1 = 1 > 0$ and ends with $y_2 = 2 > 0$ must remain in the region $y_1 > 0$ and $y_2 > 0$.

It then follows that $y_1'$ will be strictly negative, so $y_1$ is strictly decreasing with $x$, whilst $y_2'$ is strictly positive, so $y_2$ is strictly increasing with $x$. This means that you can confine your attention to $0 < y_2(0) < 2$.

 

[electronic engineer]

Are you familiar with Euler's method? It is well explained there, but I'll summarize. You know y1(0) and y2(0). From your Diff. Eq, you can calculate y1'(0) and y2'(0). Since y1' = Δy/Δx, you can write Δy1 = y1'(0) * Δx, Δy2 = y2'(0) * Δx. Then y1(Δx) = y1(0) + Δy1, y2(Δx) = y2(0) + Δy2. You just keep doing this, incrementing x by small steps (which I called Δx, and the Wiki article calls h), until you arrive at x=1.

Edit: Corrected the expressions for Δy1 and Δy2.

So when i applied your equations I got y1(1)=0 and y2(1)=2 assuming y2(0)=1 for increments h=0.2. Which is quite logical because y1 is strictly decreasing with x and y2 is strictly increasing with x as pasmith said.

 

[Chestermiller]

So when i applied your equations I got y1(1)=0 and y2(1)=2 assuming y2(0)=1 for increments h=0.2. Which is quite logical because y1 is strictly decreasing with x and y2 is strictly increasing with x as pasmith said.

Shouldn't the product of y1 and y2 be constant?

Chet

 

[electronic engineer]

Shouldn't the product of y1 and y2 be constant?

Chet

I don't think so but why do you ask?

 

[Chestermiller]

The analytic solution says that the product is constant. I ask because you got y1(1) = 0. It doesn't look like you updated the derivatives during the Euler integration.

Chet

 

[phyzguy]

So when i applied your equations I got y1(1)=0 and y2(1)=2 assuming y2(0)=1 for increments h=0.2. Which is quite logical because y1 is strictly decreasing with x and y2 is strictly increasing with x as pasmith said.

This is incorrect. Your increments are probably too large. I used h = .001. Chestermiller is right - the product y1*y2 is constant, which you can see by adding y2 *(first equation) + y1 * (second equation).

 

[electronic engineer]

This is incorrect. Your increments are probably too large. I used h = .001. Chestermiller is right - the product y1*y2 is constant, which you can see by adding y2 *(first equation) + y1 * (second equathion).

By calculating manually with the same steps for h=0.2 (which is quite large as you said) I got y2(1)=2,..(something) I don't remember but i think it's still acceptable for this large assumtption of h, so the guess that y2(0)=1 was not bad I think.

Yes and y1(1) is almost 0.5 which ensures that y1(0)*y2(0)=y1(1)*y2(1)=C=1

 

[Chestermiller]

By calculating manually with the same steps for h=0.2 (which is quite large as you said) I got y2(1)=2,..(something) I don't remember but i think it's still acceptable for this large assumtption of h, so the guess that y2(0)=1 was not bad I think.

Yes and y1(1) is almost 0.5 which ensures that y1(0)*y2(0)=y1(1)*y2(1)=C=1

The analytic solution for the initial condition on y2 you used is y2(1)= e = 2.72

Chet