visual basic excel coding (programming)

dbsq88

old_report_hw3.docx

Home >Engineering homework help >Mechanical Engineering homework help >visual basic excel coding (programming)

CALIFORNIA STATE UNIVERSITY, NORTHRIDGE

MECHANICAL ENGINEERING DEPARTMENT

MARCH 30 2015

ME 309

HOMEWORK #3

Ahmed Mohammed

Problem statement

1. Write a general computer code to solve system of up to five couples first order initial value problems using Heun and Newton iteration trapezoidal methods. These are combined in the same algor4thim in such way that Heun’s method is automatically employed to provide the initial guess at each time step for the Newton iteration required by fully implicit Heun’s method, and no iterations are needed. The complete pseudo-language algorithm is attached.

Test this code by solving the following problem.

$C:\Users\asm85688\Downloads\IMG_8677.JPG$

Use step size h= 0.1, 0.05 and 0.025. Solve the problem first with explicit Hun’s method, and then with implicit Newton iteration integration for each value of h. employ a convergence tolerance E= 0.000001 for the Newton iteration of the trapezoidal method. The exact solution to this problem is,

$C:\Users\asm85688\Downloads\IMG_8678.JPG$

Make a table of results of convergence tests (based on the exact solution) at t=1, 2, 3, 4 & 6. This table should include the exact value, the computed solution, the error and the error ratios from successive step size for each of the required values of t. discuss how these results compare with theory. Also discuss what factors should influence the choice of convergence tolerance E for the Newton iterations, and whether the specified value given above is appropriate.

2. $C:\Users\asm85688\Downloads\IMG_8679.JPG$ Solve the following problem using only the newton iteration trapezoidal method.

Employ step size h=0.1, 0.05, 0.01 and iteration convergence tolerance E= 0.000001. Consider the h=0.01 solution to be the “exact” in order to carry out convergence testes between the h= 0.1 and h=0.05 solutions at t=0.5, 1, 1.5, 2. Make a table, similar to the table in problem number 1.

Mathematical Description

HEUN’S METHOD

Heun’s Method, explained in a short manner, uses the line tangent to the function at the beginning of an interval. Now if a small step is applied to it, the error with the function result will be small. Heun’s method can be explained in more detailed in the following way:

2.-

To obtain solution point (t1,y1) we can use the fundamental theorem of calculus and integrate y’(t) over [t0,t1] to get

3.- Solving for y(t1) we find,

4.- We can use a numerical integration to approximate definite integral. If we use trapezoidal rule with step size h = t1 – t0, then we get

5.- We still need to find, y(t1) but an estimation for this value will work. After this we get the following, which is the Heun’s method.

6.- When this process is repeated it generates a sequence of points that approximate the solution curve y =y(t). At each step, Euler’s method is used as a prediction then the trapezoidal rule helps to make the correction to obtain the final value. [1]

Newton’s Method

It is way to approximate the roots of an equation by taking out the curve in the equation and then replace with a tangent line. Then, it find the intersection with the iteration process.

Suppose we want to approximate the solution to f(x) = 0 and we know the initial approximation which we call x0. To get a better approximation first we get the tangent line to f(x) at x0. So,

y = f(x0) + f’(x0)(x-x0) (from taylor series)

Now consider the following graph,

The line marking x0 is tangent to point x0. We can see that this tangent line will cross the x-axis much closer to actual solution than x0 is. At this point, we call x0 our initial guess and x1 our first approximation. To find this point x1, well, we know its coordinates (x1,0) and we know its on the tangent line, so we plug this coordinate point into the tangent formula and solve for x1.

y = f(x0) + f’(x0)(x-x0)

0= f(x0) + f’(x0)(x1-x0)

x1 = x0 –( f(x0)/f’(x0))

So we can find the new approximation provided the derivative is not zero at the original approximation. Now to get even better approximations we repeat the whole process and we the new point, x2.

X2 = x1–( f(x1)/f’(x1))

Point x2 is shown in the graph as well, it is closer to the actual answer. With the Newton’s method and a VBA code, we can approximate the answer with the accuracy we wish to obtain

Numerical Algorithm

In this assignment we asked to do 11 steps. After, we determined all the values we are going to use in the code. Starting loading the input with initial values for each question. For example, in question number 2 we used the n step equal 3, for h we use different value each time. Another example, in question number one when we use “max It” as 1 the result will be Heun’s method. If we use 10, it will be Newton iteration method. Using initial condition it will print the values on the spreadsheet. Next, we are going to begin time stepping by updating the value for t, old t and old u. This going to be the beginning of developing Newton iteration. If the program found the value of max it “number of iteration” bigger than 1 it will step to the Jacobin function. If max it is equal to on it will continue to the next step and then it will jump to step 11 where the program will print the result. For question number one the program will go to step 6 and it will print the result in the spreadsheet as Heun’s method.

However, in question number 2 it is going to be different. The program will complete Heun’s method for use as initial guess for trapezoidal rule. It will solve for trapezoidal rule. It will calculate it using forward Euler, where it’s going to obtaining an approximate solution of an ordinary differential equation. Evaluating the area between a curve and axis by approximate the area with the area of trapezoidal. Then, it will load Jacobian for Newton’s iteration. The program will also going to update the coefficient matrix with the Jacobian. Using Gaussian elimination it will solve for delta U and test convergence of Newton iteration. Lastly, it will print out the result on the spreadsheet with the values of “U”, time and number of iteration.

In the code we had to do some function sub so the program can call it when its need it. We made derivative function for each question, Jacobin function and Gaussian Elimination.

Production run and Success criteria

Problem #1

$U:\ME309\Capture1.PNG$

$U:\ME309\Capture2.PNG$

As you can see form the result above we can compare the values of exact value, calculated value and the error ratio. Heun’s method looks closer to the actual result of Uexact and we can know that from the ratio where it’s much smaller than newton value.

Problem #2

$U:\ME309\Capturelast.PNG$

Calculated value when it is smaller it get closer to the value as you can see the difference between h=0.05 and 0.1 to the Uexact value. At h=0.05 is more accurate. However, the goal from calculating the error ratio is to see how far we got from the exact value at certain point. Also, we stopped approximating the answer when we were as close as 1x10-6.

Program Code

Problem #2

Sub NEWTON()

'1 load input data

Dim neqns As Double

Dim maxit As Double

Dim nsteps As Double

Dim u(5) As Double

Dim t, n, i As Double

Dim f(5) As Double

Dim h As Double

Dim uold(5) As Double

Dim fold(5) As Double

Dim ustar(5) As Double

Dim A(5, 5) As Double

Dim fupp(5) As Double

Dim delta(5, 5) As Double

Dim deltau(5) As Double

Dim deltaumax As Double

u(1) = 0

u(2) = 0

u(3) = 0

E = 0.000001

i = 1

h = 0.05

neqns = 3

maxit = 10

nsteps = 2 / h

'2: begin time stepping and update

For n = 1 To nsteps

told = t

t = told + h

For i = 1 To neqns

uold(i) = u(i)

Next i

'3: begin newton iterations

For m = 1 To maxit

If m > 1 Then GoTo 6

'4: evaluate u*(i) for hun's method

For i = 1 To neqns

fold(i) = fderv(uold, told, i)

ustar(i) = uold(i) + h * fold(i)

Next i

amax = Abs(A(i, k))

rowmax = i

End If

Next i

'swap Rows

For j = k To n + 1

temp = A(k, j)

A(k, j) = A(rowmax, j)

A(rowmax, j) = temp

Next j

'MAIN GAUSS FUNCTION

For i = k + 1 To n

p = A(i, k) / A(k, k)

For j = k To n + 1

A(i, j) = A(i, j) - A(k, j) * p

Next j

Next i

Next k

'OUTPUT

x(n) = A(n, n + 1) / A(n, n)

For i = n - 1 To 1 Step -1

SumProduct = 0

For j = i + 1 To n

SumProduct = A(i, j) * x(j) + SumProduct

Next j

x(i) = (A(i, n + 1) - SumProduct) / A(i, i)

Next i

End Sub

5: ' complate heun's method foruse as initial guess for trapezoidal rule

For i = 1 To neqns

f(i) = fderv(ustar, t, i)

u(i) = uold(i) + 0.5 * h * (fold(i) + f(i))

Next i

If maxit = 1 Then GoTo 11

6: ' load j(f) into the cofficient matrix (A)

Call jacobian(A, neqns, u, t)

7: 'Evaluate F(u) [for the last column of the Augm Coeff

For i = 1 To neqns

f(i) = fderv(u, t, i)

fupp(i) = u(i) - 0.5 * h * f(i) - (uold(i) + 0.5 * h * fold(i))

A(i, neqns + 1) = -fupp(i)

For j = 1 To neqns

If i = j Then delta(i, j) = 1

If i <> j Then delta(i, j) = 0

A(i, j) = delta(i, j) - (0.5 * h * A(i, j))

Next j

Next i

8: 'Solve for deltau(i) using gaussian elimination

Call gauss(A, neqns, deltau)

9: 'Calculate abs(deltau) and increment u(i),

deltaumax = 0

For i = 1 To neqns

If Abs(deltau(i)) > deltaumax Then deltaumax = Abs(deltau(i))

u(i) = u(i) + deltau(i)

Next i

10: ' test convergence of newton iteration

If deltaumax < E Then GoTo 11

Next m

11 ': print result for u(i)

Cells(n + 1, 1).Value = t

Cells(n + 1, 2).Value = u(1)

Cells(n + 1, 3).Value = u(2)

Cells(n + 1, 4).Value = u(3)

Cells(n + 1, 5).Value = m

Next n

End Sub

Function fderv(u, t, i)

If i = 1 Then fderv = u(2)

If i = 2 Then fderv = u(3)

If i = 3 Then fderv = -2 * u(1) * u(3) + ((u(2) ^ 2) - 3 * u(1) + 10 * Sin(6 * t))

End Function

Function jacobian(A, neqns, u, t)

A(1, 1) = 0

A(1, 2) = 1

A(1, 3) = 0

A(2, 1) = 0

A(2, 2) = 0

A(2, 3) = 1

A(3, 1) = -3 - 2 * u(3)

A(3, 2) = 2 * u(2)

A(3, 3) = -2 * u(1)

End Function

Sub gauss(A, n, x)

For k = 1 To n - 1

'COLUMN PIVOTING

'identify largest abs value in cell

amax = Abs(A(k, k))

rowmax = k

For i = k + 1 To n

If Abs(A(i, k)) > amax Then

Discussion

In this assignment we learned how to solve differential equations using two new method. I notice that Heun’s method came out to be more accurate than Newton method. From calculating the ratio for each method as shown in the tables we conclude Heun’s method is more accurate. However, when comparing the exact value we found the values for Heun’s method is closer. Newton method dose the work for approximate solutions to equations as Heun’s method but it is not that accurate. The basic idea behind these processes is to take advantage of what a tangent line is in the form of an equation, and use this same equation with the new values until the answer is close enough to the actual answer.

Reference

[1] http://mathfaculty.fullerton.edu/mathews//n2003/heunsmethod/Heun%27sMethodProof.pdf

[2] http://mathworld.wolfram.com/NewtonsMethod.html