a8daaf7b-1b1d-4f46-8011-980372ff544d 33 15 12 decimal

&[DATE] &[TIME] - &[FILENAME]

Adapted Mathcad: Nathan Collier, Autar Kaw University of South Florida

Section 1: input

Differential Equation to illustrate the method A1 y x 2 diff * A2 y x diff * + A3 y * + A4 + 0 ≡

Boundary Cnditions x0 0 : y0 10 : xf 1 : yf 0 : 41 sys

Explicit/Value the system x y A1 1 : x y A2 1 : x y A3 1 : x A4 1 : 41 sys

initial values dydx1 15 - : dydx2 10 - : 21 sys

Section 2: Setup of the Approximate Solution

To use RK4method of solving ODEs, the equations must be of order one. Thus, the derivative dy/dx is set to another variable and a system of first order is made.

y z ≡ x diff

z x diff A2 z * A3 y * + A4 + A1 - / ≡

These two equations are now in the proper form
to apply RK4 method. We need to define two 
functions to be used in determining 'k' values. x y z f1 z : x y z f2 x y A2 z * x y A3 y * + x A4 + x y A1 - / : 21 sys

The step size will be used in this document h xf x0 - 5 / : 0.2

Section 3: First Iteration [the 5 steps]

Step 1

X1 1 el x0 : Y1 1 el y0 : Z1 1 el dydx1 : X1 2 el X1 1 el h + : 41 line

X1 Y1 Z1 31 mat 0 0.2 21 mat 1011 mat 15 - 11 mat 31 mat

k1y X1 1 el Y1 1 el Z1 1 el f1 : 15 -

k1z X1 1 el Y1 1 el Z1 1 el f2 : 4

k2y X1 1 el 0.5 h * + Y1 1 el 0.5 h * k1y * + Z1 1 el 0.5 h * k1z * + f1 : 14.6 -

k2z X1 1 el 0.5 h * + Y1 1 el 0.5 h * k1y * + Z1 1 el 0.5 h * k1z * + f2 : 5.1

k3y X1 1 el 0.5 h * + Y1 1 el 0.5 h * k2y * + Z1 1 el 0.5 h * k2z * + f1 : 14.49 -

k3z X1 1 el 0.5 h * + Y1 1 el 0.5 h * k2y * + Z1 1 el 0.5 h * k2z * + f2 : 4.95

k4y X1 1 el h + Y1 1 el h k3y * + Z1 1 el h k3z * + f1 : 14.01 -

k4z X1 1 el h + Y1 1 el h k3y * + Z1 1 el h k3z * + f2 : 5.908

Y1 2 el Y1 1 el h 6 / k1y 2 k2y * + 2 k3y * + k4y + (* + : 10 7.094 21 mat

Z1 2 el Z1 1 el h 6 / k1z 2 k2z * + 2 k3z * + k4z + (* + : 15 - 14 - 21 mat

Step 2

X1 3 el X1 2 el h + : 0 0.2 0.4 31 mat

k1y X1 2 el Y1 2 el Z1 2 el f1 : 14 -

k1z X1 2 el Y1 2 el Z1 2 el f2 : 5.906

k2y X1 2 el 0.5 h * + Y1 2 el 0.5 h * k1y * + Z1 2 el 0.5 h * k1z * + f1 : 13.409 -

k2z X1 2 el 0.5 h * + Y1 2 el 0.5 h * k1y * + Z1 2 el 0.5 h * k1z * + f2 : 6.715

k3y X1 2 el 0.5 h * + Y1 2 el 0.5 h * k2y * + Z1 2 el 0.5 h * k2z * + f1 : 13.328 -

k3z X1 2 el 0.5 h * + Y1 2 el 0.5 h * k2y * + Z1 2 el 0.5 h * k2z * + f2 : 6.575

k4y X1 2 el h + Y1 2 el h k3y * + Z1 2 el h k3z * + f1 : 12.685 -

k4z X1 2 el h + Y1 2 el h k3y * + Z1 2 el h k3z * + f2 : 7.257

Y1 3 el Y1 2 el h 6 / k1y 2 k2y * + 2 k3y * + k4y + (* + : 10 7.094 4.422 31 mat

Z1 3 el Z1 2 el h 6 / k1z 2 k2z * + 2 k3z * + k4z + (* + : 15 - 14 - 12.675 - 31 mat

Step 3

X1 4 el X1 3 el h + : 0 0.2 0.4 0.6 41 mat

k1y X1 3 el Y1 3 el Z1 3 el f1 : 12.675 -

k1z X1 3 el Y1 3 el Z1 3 el f2 : 7.253

k2y X1 3 el 0.5 h * + Y1 3 el 0.5 h * k1y * + Z1 3 el 0.5 h * k1z * + f1 : 11.95 -

k2z X1 3 el 0.5 h * + Y1 3 el 0.5 h * k1y * + Z1 3 el 0.5 h * k1z * + f2 : 7.795

k3y X1 3 el 0.5 h * + Y1 3 el 0.5 h * k2y * + Z1 3 el 0.5 h * k2z * + f1 : 11.895 -

k3z X1 3 el 0.5 h * + Y1 3 el 0.5 h * k2y * + Z1 3 el 0.5 h * k2z * + f2 : 7.669

k4y X1 3 el h + Y1 3 el h k3y * + Z1 3 el h k3z * + f1 : 11.141 -

k4z X1 3 el h + Y1 3 el h k3y * + Z1 3 el h k3z * + f2 : 8.099

Y1 4 el Y1 3 el h 6 / k1y 2 k2y * + 2 k3y * + k4y + (* + : 10 7.094 4.422 2.038 41 mat

Z1 4 el Z1 3 el h 6 / k1z 2 k2z * + 2 k3z * + k4z + (* + : 15 - 14 - 12.675 - 11.132 - 41 mat

Step 4

X1 5 el X1 4 el h + :

k1y X1 4 el Y1 4 el Z1 4 el f1 : 11.132 -

k1z X1 4 el Y1 4 el Z1 4 el f2 : 8.094

k2y X1 4 el 0.5 h * + Y1 4 el 0.5 h * k1y * + Z1 4 el 0.5 h * k1z * + f1 : 10.323 -

k2z X1 4 el 0.5 h * + Y1 4 el 0.5 h * k1y * + Z1 4 el 0.5 h * k1z * + f2 : 8.398

k3y X1 4 el 0.5 h * + Y1 4 el 0.5 h * k2y * + Z1 4 el 0.5 h * k2z * + f1 : 10.292 -

k3z X1 4 el 0.5 h * + Y1 4 el 0.5 h * k2y * + Z1 4 el 0.5 h * k2z * + f2 : 8.287

k4y X1 4 el h + Y1 4 el h k3y * + Z1 4 el h k3z * + f1 : 9.475 -

k4z X1 4 el h + Y1 4 el h k3y * + Z1 4 el h k3z * + f2 : 8.495

Y1 5 el Y1 4 el h 6 / k1y 2 k2y * + 2 k3y * + k4y + (* + : 10 7.094 4.422 2.038 0.023 - 51 mat

Z1 5 el Z1 4 el h 6 / k1z 2 k2z * + 2 k3z * + k4z + (* + : 15 - 14 - 12.675 - 11.132 - 9.467 - 51 mat

Step 5

X1 6 el X1 5 el h + :

X1 0 0.2 0.4 0.6 0.8 161 mat

k1y X1 5 el Y1 5 el Z1 5 el f1 : 9.467 -

k1z X1 5 el Y1 5 el Z1 5 el f2 : 8.49

k2y X1 5 el 0.5 h * + Y1 5 el 0.5 h * k1y * + Z1 5 el 0.5 h * k1z * + f1 : 8.618 -

k2z X1 5 el 0.5 h * + Y1 5 el 0.5 h * k1y * + Z1 5 el 0.5 h * k1z * + f2 : 8.588

k3y X1 5 el 0.5 h * + Y1 5 el 0.5 h * k2y * + Z1 5 el 0.5 h * k2z * + f1 : 8.608 -

k3z X1 5 el 0.5 h * + Y1 5 el 0.5 h * k2y * + Z1 5 el 0.5 h * k2z * + f2 : 8.493

k4y X1 5 el h + Y1 5 el h k3y * + Z1 5 el h k3z * + f1 : 7.768 -

k4z X1 5 el h + Y1 5 el h k3y * + Z1 5 el h k3z * + f2 : 8.513

Y1 6 el Y1 5 el h 6 / k1y 2 k2y * + 2 k3y * + k4y + (* + : 10 7.094 4.422 2.038 0.023 - 1.746 - 61 mat

Z1 6 el Z1 5 el h 6 / k1z 2 k2z * + 2 k3z * + k4z + (* + : 15 - 14 - 12.675 - 11.132 - 9.467 - 7.761 - 61 mat

x0 0 : y0 10 : xf 1 : yf 0 : 41 sys

BoundaryConditions 01012 mat 1012 mat 21 mat :

iterX X1 Y1 augment :

iterX 010 0.2 7.094 0.4 4.422 0.6 2.038 0.8 0.023 - 1 1.746 - 62 mat

data char size clr plot k 1 data rows range r3 k el char : r4 k el size : r5 k el clr : 13 mat for data r3 r4 r5 augment 21 line :

BC_start BoundaryConditions 1 el . 14 black plot : BC_end BoundaryConditions 2 el . 14 black plot : iter5 iterX o 7 red plot : 31 sys

BC_start BC_end iter5 3 1 sys

Section 4: Second Iteration [forloop]

is defined for all from the last 5th iteration X X1 :

Y2 Z 1 el dydx2 : Y 1 el y0 : i 15 range k1y X i el Y i el Z i el f1 : k1z X i el Y i el Z i el f2 : k2y X i el 0.5 h * + Y i el 0.5 h * k1y * + Z i el 0.5 h * k1z * + f1 : k2z X i el 0.5 h * + Y i el 0.5 h * k1y * + Z i el 0.5 h * k1z * + f2 : k3y X i el 0.5 h * + Y i el 0.5 h * k2y * + Z i el 0.5 h * k2z * + f1 : k3z X i el 0.5 h * + Y i el 0.5 h * k2y * + Z i el 0.5 h * k2z * + f2 : k4y X i el h + Y i el h k3y * + Z i el h k3z * + f1 : k4z X i el h + Y i el h k3y * + Z i el h k3z * + f2 : Y i 1 + el Y i el h 6 / k1y 2 k2y * + 2 k3y * + k4y + (* + : Z i 1 + el Z i el h 6 / k1z 2 k2z * + 2 k3z * + k4z + (* + : 101 line for Y 41 line :

iter2 X Y2 augment :

iter2 010 0.2 7.994 0.4 6.027 0.6 4.162 0.8 2.449 1 0.922 62 mat

BC_start BoundaryConditions 1 el . 14 black plot : BC_end BoundaryConditions 2 el . 14 black plot : iter5 iterX o 7 red plot : ITER2 iter2 o 7 blue plot : 41 sys

BC_start BC_end iter5 ITER2 4 1 sys

Section 5: Third Iteration [forloop]

By examining the figure immediately above,you will note that at 'x=x0'the solution [for both iterations] is the boundary condition. At 'x=xf' however, it is likely that the solution is NOT, unless you are a good guesser. We can use the first two iterations to predict what a better value for dy/dx should be. We will assume that the relationship is linear and use interpolation to determine the final estimate for dy/dx. This will yield exact results if the equation is linear while for non-linear differential equations it will be a reasonable approximation.

The third assumed dy/dx will be dydx3 dydx1 dydx2 - (Y1 6 el Y2 6 el - (/ yf Y2 6 el - (* dydx2 + : 11.727 -

Y3 Z 1 el dydx3 : Y 1 el y0 : i 15 range k1y X i el Y i el Z i el f1 : k1z X i el Y i el Z i el f2 : k2y X i el 0.5 h * + Y i el 0.5 h * k1y * + Z i el 0.5 h * k1z * + f1 : k2z X i el 0.5 h * + Y i el 0.5 h * k1y * + Z i el 0.5 h * k1z * + f2 : k3y X i el 0.5 h * + Y i el 0.5 h * k2y * + Z i el 0.5 h * k2z * + f1 : k3z X i el 0.5 h * + Y i el 0.5 h * k2y * + Z i el 0.5 h * k2z * + f2 : k4y X i el h + Y i el h k3y * + Z i el h k3z * + f1 : k4z X i el h + Y i el h k3y * + Z i el h k3z * + f2 : Y i 1 + el Y i el h 6 / k1y 2 k2y * + 2 k3y * + k4y + (* + : Z i 1 + el Z i el h 6 / k1z 2 k2z * + 2 k3z * + k4z + (* + : 101 line for Y 41 line :

iter3 X Y3 augment :

iter3 010 0.2 7.683 0.4 5.472 0.6 3.428 0.8 1.595 1062 mat

BC_start BoundaryConditions 1 el . 20 black plot : BC_end BoundaryConditions 2 el . 20 black plot : iter5 iterX . 14 red plot : ITER2 iter2 . 14 blue plot : ITER3 iter3 . 14 green plot : 51 sys

BC_start BC_end iter5 ITER2 ITER3 5 1 sys

iterX iter2 iter3 13 mat 010 0.2 7.094 0.4 4.422 0.6 2.038 0.8 0.023 - 1 1.746 - 62 mat 010 0.2 7.994 0.4 6.027 0.6 4.162 0.8 2.449 1 0.922 62 mat 010 0.2 7.683 0.4 5.472 0.6 3.428 0.8 1.595 1062 mat 13 mat