Numerical Analysis and Simulation

Since the questioner didn't specify which shooting method was being used, I'll explain the linear shooting method.

First off, the shooting method is used to solve boundary-value problems, or BVPs, where instead of being given 2 initial conditions to solve a problem, we only have the value of the function at both ends of a given interval [a,b], i.e. we are given y(a) and y(b) as conditions for a problem y''=f(x,y,y'), where x is in the interval [a,b].

Now, suppose our problem is written in the form of y''= p(x)y'+q(x)y+r(x), y(a)=s and y(b)=t with the functions p, q and r continuous on the interval [a,b] and q(x)>0 on [a,b], I wont go into the theory of it but this means we can separate our problem into a system of 2 initial-value problems, namely:

y1''=p(x)y1'+q(x)y1+r(x), x in [a,b], y(a)=s, y'(a)=0

y2''=p(x)y2'+q(x)y2, x in [a,b], y(a)=0, y'(a)=1

and that the linear combination y(x)=y1(x)+((t-y1(b))/y2(b))*y2(a) is a unique solution to our linear boundary-value problem.

The textbook I am using for the majority of this information gives the following algorithm for the shooting method with Runge-Kutta 4 in terms of its individual components, in order to save a LOT of space, I'm gonna try and write it up in terms of vectors instead.

The algorithm for the shooting method with R-K4 approximating y1 and y2 is as follows:

INPUT: endpoints a, b; boundary conditions s, t; number of sub intervals N; functions p, q and r of x; vectors u=(s, 0) v=(0,1) (please note that in the algorithm a vector followed by [1] or [2] signifies the first/second component of the vector. e.g. at the moment u[1]=s and u[2]=0)

STEP 1: set h=(b-a)/N

STEP 2: For i=0...N-1 do STEPS 3 to 5

STEP 3: set x=a+ih

STEP 4: set;

k1=(hu[2] , h[p(x)u[2]+q(x)u[1]+r(x)])

k2=(hu[2]+0.5*k1[2] , h[p(x+h/2)*(u[2]+0.5*k1[2])+q(x+h/2)*(u[1]+0.5*k1[1])+r(x+h/2)])

k3=(hu[2]+0.5*k2[2] , h[p(x+h/2)*(u[2]+0.5*k2[2])+q(x+h/2)*(u[1]+0.5*k2[1])+r(x+h/2)])

k4=(hu[2]+k3[2] , h[p(x+h)*(u[2]+k3[2])+q(x+h)*(u[1]+k3[1])+r(x+h)])

set u[1]=(u[1]+(k1[1]+2*k2[1]+2*k3[1]+k4[1])/6)

set u[2]=(u[2]+(k1[2]+2*k2[2]+2*k3[2]+k4[2])/6)

STEP 5: set;

k'1=(hv[2] , h[p(x)v[2]+q(x)v[1]])

k'2=(hv[2]+0.5*k'1[2] , h[p(x+h/2)*(v[2]+0.5*k'1[2])+q(x+h/2)*(v[1]+0.5*k'1[1])])

k'3=(hv[2]+0.5*k'2[2] , h[p(x+h/2)*(v[2]+0.5*k'2[2])+q(x+h/2)*(v[1]+0.5*k'2[1])])

k'4=(hv[2]+k'3[2] , h[p(x+h)*(v[2]+k'3[2])+q(x+h)*(v[1]+k'3[1])])

set v[1]=(v[1]+(k'1[1]+2*k'2[1]+2*k'3[1]+k'4[1])/6)

set v[2]=(v[2]+(k'1[2]+2*k'2[2]+2*k'3[2]+k'4[2])/6)

STEP 6: set w=(s, (t-u[1])/v[1])

OUTPUT (a, w)

STEP 7: for i=1...N;

set W1=u[1,i]+w[2]*v[1,i]

set W2=u[2,i]+w[2]*v[2,i]

set x+a+ih;

OUTPUT (x, W1, W2) (output is actually components of the ith iteration of the vector w, it's not giving out the last iteration but a list of all the outputs as the last iteration isn't very helpful for a boundary problem, we want to know the values inbetween)

STEP 8: STOP

The value of the function y(x) is then given by y(xi)=y1(xi)+y2(xi)*(t-y1(b))/y2(b)

Hopefully you can read the algorithm, just for clarification on the running of it, it does one run on i=0 for the R-K method (steps 4 and 5), then again for i=1, it then moves to the shooting method (step 7) and does one iteration using the vector values generated by the R-K4 on that run to generate the value of y(x) on the first step after x=a i.e. at x=a+h. It then goes to i=2 and does a run of R-K generating new vector values again, and then goes to the shooting method and generates the next value of the function at x=a+2h and so on until it has finished generating a list of values for W1 and W2 which correspond to y1(xi) and y2(xi) respectively. You can then use this values with the formula above to find the values of the function at the step points between a and b.

However in actuality the algorithm doesn't run like that, I just think it easier to think of it in that way. What it actually does is run step 3-5 for i=0 to i=N-1, so on the final iteration you generate a value for u and v at i=N. you use these Nth iterations in step 6. step 7 then runs from i=1 to i=N using the vector values previously generated for the corresponding value of i. Hence why I felt then need to put u[1,i] etc. in, in hopes of clarifying what value the algorithm is using.

And there you have the shooting method with a Runge-Kutta 4 technique for linear BVPs. Hopefully this made sense, if not send me a message detailing what it is you don't understand and I will see if I can help if I have spare time.

This is gonna turn into a long answer, please bear with me. If you already know the algorithm's just skip to the examples.

Iterative techniques for solving a linear system Ax=b,where A is a square matrix with non-zero diagonal entries, start of by first rearranging the system into a form x=Tx+c where T is a fixed square matrix with zero entries along the diagonal. an initial approximation x(0) is then applied to the Tx+c part and gives the next set of of approximation values for x(1). First I'll give an example of the Jaboci method and then the Gauss-Seidal method. Both examples can be found at the start of section 7.3 of Numerical Analysis (8th edition) by R. L. Burden and J. Douglas Faires.

The Jacobi method generates the next approximation by using the formula,

xi(k+1)=Σnj=1, j=/=i (-aijxj(k))/aii + bi/aii.

Note that aii must not be equal to zero which is why we have A must be a square matrix with non-zero diagonal entries.

To show what I mean we will use the following system:

E1= 10x1 - x2 + 2x3 = 6

E2= -x1 + 11x2 - x3 + 3x4 =25

E3= 2x1 - x2 + 10x3 -x4 = -11

E4= 3x2- x3 + 8x4 = 15

Converting this from the form Ax=b to x=Tx+c we get:

x1(k+1) = x2(k)/10 - x3(k)/5 +3/5

x2(k+1) = x1(k)/11 + x3(k)/11 - 3x4(k)/11 + 25/11

x3(k+1) = -x1(k)/5 + x2(k)/10 + x4(k)/10 - 11/10

x4(k+1) = -3x2(k)/8 + x3(k)/8 + 15/8

Notice how this forms a matrix T where the diagonals are all zero and c=(3/5, 25/11, -11/10, 15/8).

For an initial approximation we shall use x(0)=(0,0,0,0) plugging these values we into the Tx+c part we get x(1)=(0.6, 2.2727, -1.1, 1.8750) and plugging these values in we get x(2)=(1.0473, 1.7159, -0.8052, 0.8852). After 10 iterations we obtain x(10)=(1.0001, 1.9998, -0.9998, 0.9998) giving us a maximum absolute error of 0.0002 as the unique solution is x=(1, 2, -1, 1)

You may notice that it does not matter in which order you compute the equations as all the variable only take the previous iteration. This is where The Gauss-Seidal method improves upon the Jacobi method to make a 'better' iteration method.

To explain it better, go back to the iteration equations above, if you were to do them in order then for x2(k+1)=... part the first term you have is x1(k)/11, but in the previous equation you just worked out x(k+1). Surely it makes more sense to use this new value in the iteration for x2(k+1) than the old value, and this is exactly what the G-S method does! The G-S method generates the next approximation using the formula:

xi(k+1) = (-Σi-1j=1(aijxj(k+1)) - Σnj=i+1(aijxj(k)) + bi)/aii

Using our previous example this changes our iteration equations into the following:

x1(k+1) = x2(k)/10 - x3(k)/5 +3/5

x2(k+1) = x1(k+1)/11 + x3(k)/11 - 3x4(k)/11 + 25/11

x3(k+1) = -x1(k+1)/5 + x2(k+1)/10 + x4(k)/10 - 11/10

x4(k+1) = -3x2(k+1)/8 + x3(k+1)/8 + 15/8

Again starting with an initial approximation x(0)=(0, 0, 0, 0) and plugging into the equations above we get x(1) = (0.6, 2.3272, -0.9873, 0.8789) and now put these into the equations above and repeat. This time after 5 iterations we get x(5)=(1.0001, 2, -1, 1) so the G-S method has given us an answer with max absolute error of 0.0001 in half the amount of iterations that it took the Jacobi method to get an answer with max absolute error of 0.0002!

The answer depends on what information you have. Given only one number, all that can be said is that the base is larger than the largest digit appearing in the number.

Equations will not help if there are no "carries". For example,

For example, 10 + 12 = 22 is true in any base greater than or equal to 3.

[In base 3, it is 3+5=8, in base 4, it is 4+6=10, in base 500 it is 500+502=1002, in base e (= 2.7183...) it is 2.7183+4.7183 =7.4366 (approx).]

Given an infinitely convergent sequence pn with limit p, the forward difference is the measure of the difference between the current term and he next. The backward difference is the measure of the difference between the current term and the previous.

i.e.

forward difference: Δpn=pn+1 - pn

backward difference: ∇pn=pn - pn-1

Also, note that since they are both expressed by pn, the forward difference is recognised by the use of a delta before the pn, and the backward difference by the use of a nabla.