[Solved] AMATH301 Homework 3

1. Consider the circuit diagram below:

Here, each V represents a change in voltage (measured in volts) at a battery, each R represents a resistance (measured in ohms) at a resistor, and each I represents a current (measured in amps) through a wire. These quantities obey two simple

laws:

• Ohm’s Law: The voltage drop across a resistor is V = IR,
• Kirchhoff’s Second Law: The sum of all the voltage drops in a closed loop

is zero.

Using these two laws, we can construct the following systems of equations:

R₆I₁ + R₁ (I₁ − I₂) + R₂ (I₁ − I₃) =V₁,

R₃I₂ + R₄ (I₂ − I₃) + R₁ (I₂ − I₁) =V₂, R₅I₃ + R₄ (I₃ − I₂) + R₂ (I₃ − I₁) =V₃.

Let the resistances be given by R₁ = 10Ω, R₂ = 20Ω, R₃ = 5Ω, R₄ = 15Ω, R₅ = 30Ω and R₆ = 25Ω. We want to calculate the currents I₁, I₂, and I₃.

• Write the equations in matrix form Ax = b, where x is a 3 × 1 vector of the unknown current values (you need to do this by hand, not in Matlab). Using the lu command in Matlab, find matrices L, U and P such that PA = LU. Calculate the product UPL and save the resulting 3 × 3 matrix in dat. (Notice that I haven’t told you the voltages yet. Does that matter?)
• Let V₁ = 50 and V₂ = 0. For every value of V₃ from 1 to 100 in increments of

1, calculate I₁, I₂ and I₃ using LU decomposition. Save the resulting values of I₂ as a 1 × 100 row vector in A2.dat.

• Repeat part (b), but solve each system using the inv command (i.e. using the inverse of A) instead of LU decomposition. This method should be slower, and the result should be slightly different. Similar to part (b), make 1 × 100 row vector of the values of I₂. Subtract your vector from part (c) from the vector from part (b), and save the difference as dat. For example, if x_b is the name of your vector from part (b) and x_c is the name of your vector from part (c), then the vector you should save in A3.dat is x_b – x_c.

Things to think about: The inv command is not very useful in practical applications. In fact, Matlab probably warns you not to use it in your code. However, you probably saw in this code that it does not make much of a difference to the final answer. Can you find a system Ax = b where the inv command causes a lot of rounding error? Can you find a system where the inv command is visibly slower than backslash?

As we learned in class, the backslash command uses LU decomposition for most systems. Why was it better to find L, U, and P ourselves before doing part

(b)?

2. Consider the bridge truss diagrammed below.

 −s −s

10−10000000000

001 0−100000000

0000 s−s−s000000

000000−1100000

0000000−101000

00000000 0−1010

0000000000 0−10

0−100000 000001

s−s00000 0 0−ss00

000−1000 000100

000 0−s−s0 0 0 s s00

00 0 0 0  −1  0 0 1 0 0 0  0

A = 

00000000000

Given a vector of external forces b on the bridge, we can compute the forces F₁, F₂ ,…,F₁₃ by solving the system Ax = b, where x is a vector of unknown forces and A is given by

√

and s = 2/2. We will solve for the vector of forces x assuming that there are three vehicles at positions 6, 7 and 8 with weights 5 tons, 10 tons, and 4 tons, respectively. This means that b = [0,0,0,0,0,0,0,0,4,0,10,0,5]^T. (Note that the

T means transpose so b is a 13 × 1 column vector with those entries.)

• Solve for x using LU decomposition. Save the intermediate answer y in dat and the final answer x in A5.dat.
• Solve for x using the backslash command. Save your answer as dat. How does this compare with your answer from part (a)?
• Now suppose that we add weight to the truck at position 6 (which corresponds to the 13th entry of b) in increments of 0.01 tons until the bridge collapses. Each bridge member can withstand 30 tons of compression or tension (i.e., positive or negative forces.) Therefore, the bridge will collapse when the absolute value of the largest force is larger than 30. Find the smallest weight of the truck at position 6 for which the bridge collapses. Save your answer as dat. Hint:

You may find some of the functions max, abs, or norm useful.

3. Consider the matrix

.

Find the condition number of A and save it in A8.dat. You should find that the condition number is not particularly large, which suggests that any reasonable algorithm should work well with this matrix. It is easy to check by hand that an

LU decomposition of A is

Multiply LU by hand and confirm that A = LU. Now use Matlab to multiply LU and save your answer in A9.dat. Is this close to A?

If we switch the rows of A, we get a new matrix

.

An LU decomposition of B is

.

Multiply LU by hand and confirm that B = LU. Now use Matlab to multiply LU and save your answer in A10.dat. Is this close to B?

Things to think about: Why do you think one of these answers is so inaccurate? (This requires some knowledge about how floating point numbers are rounded.) What happens if you try to use the L and U matrices from above to solve Ax = b? In particular, try b = [1 + 10⁻²⁰,2]^T. The solution to that system is x = [1,1]^T.) Does the same issue arise with the matrix B? Try using the lu command to do the LU decomposition of A. Does Matlab give the same LU decomposition as above? How does Matlab’s LU decomposition of A compare with the LU decomposition of B above?