# tgr12001 / ME3255S2017

forked from sed12008/ME3255S2017
Older
100644 699 lines (699 sloc) 29.2 KB
Feb 27, 2017
1
{
2
"cells": [
3
{
4
"cell_type": "code",
5
"execution_count": 18,
6
7
"collapsed": true
8
},
9
"outputs": [],
10
"source": [
11
"%plot --format svg"
12
]
13
},
14
{
15
"cell_type": "code",
16
"execution_count": 19,
17
18
"collapsed": true
19
},
20
"outputs": [],
21
"source": [
22
"setdefaults"
23
]
24
},
25
{
26
"cell_type": "markdown",
27
28
"source": [
29
"What is the determinant of A?\n",
30
"\n",
31
"\$A=\\left[ \\begin{array}{ccc}\n",
32
"10 & 2 & 1 \\\\\n",
33
"0 & 1 & 1 \\\\\n",
34
"0 & 0 & 10\\end{array} \\right]\$\n",
35
"\n",
36
"![responses to determinant of A](det_A.png)"
37
]
38
},
39
{
40
"cell_type": "markdown",
41
42
"source": [
43
"# LU Decomposition\n",
44
"### efficient storage of matrices for solutions\n",
45
"\n",
46
"Considering the same solution set:\n",
47
"\n",
48
"\$y=Ax\$\n",
49
"\n",
50
"Assume that we can perform Gauss elimination and achieve this formula:\n",
51
"\n",
52
"\$Ux=d\$ \n",
53
"\n",
54
"Where, \$U\$ is an upper triangular matrix that we derived from Gauss elimination and \$d\$ is the set of dependent variables after Gauss elimination. \n",
55
"\n",
56
"Assume there is a lower triangular matrix, \$L\$, with ones on the diagonal and same dimensions of \$U\$ and the following is true:\n",
57
"\n",
58
"\$L(Ux-d)=Ax-y=0\$\n",
59
"\n",
60
"Now, \$Ax=LUx\$, so \$A=LU\$, and \$y=Ld\$.\n",
61
"\n",
62
"\$2x_{1}+x_{2}=1\$\n",
63
"\n",
64
"\$x_{1}+3x_{2}=1\$\n",
65
"\n",
66
"\n",
67
"\$\\left[ \\begin{array}{cc}\n",
68
"2 & 1 \\\\\n",
69
"1 & 3 \\end{array} \\right]\n",
70
"\\left[\\begin{array}{c} \n",
71
"x_{1} \\\\ \n",
72
"x_{2} \\end{array}\\right]=\n",
73
"\\left[\\begin{array}{c} \n",
74
"1 \\\\\n",
75
"1\\end{array}\\right]\$\n",
76
"\n",
77
"f21=0.5\n",
78
"\n",
79
"A(2,1)=1-1 = 0 \n",
80
"\n",
81
"A(2,2)=3-0.5=2.5\n",
82
"\n",
83
"y(2)=1-0.5=0.5\n",
84
"\n",
85
"\$L(Ux-d)=\n",
86
"\\left[ \\begin{array}{cc}\n",
87
"1 & 0 \\\\\n",
88
"0.5 & 1 \\end{array} \\right]\n",
89
"\\left(\\left[ \\begin{array}{cc}\n",
90
"2 & 1 \\\\\n",
91
"0 & 2.5 \\end{array} \\right]\n",
92
"\\left[\\begin{array}{c} \n",
93
"x_{1} \\\\ \n",
94
"x_{2} \\end{array}\\right]-\n",
95
"\\left[\\begin{array}{c} \n",
96
"1 \\\\\n",
97
"0.5\\end{array}\\right]\\right)=0\$\n"
98
]
99
},
100
{
101
"cell_type": "code",
102
"execution_count": 3,
103
104
"collapsed": false
105
},
106
"outputs": [
107
{
108
"name": "stdout",
109
"output_type": "stream",
110
"text": [
111
"A =\n",
112
"\n",
113
" 2 1\n",
114
" 1 3\n",
115
"\n",
116
"L =\n",
117
"\n",
118
" 1.00000 0.00000\n",
119
" 0.50000 1.00000\n",
120
"\n",
121
"U =\n",
122
"\n",
123
" 2.00000 1.00000\n",
124
" 0.00000 2.50000\n",
125
"\n",
126
"ans =\n",
127
"\n",
128
" 2 1\n",
129
" 1 3\n",
130
"\n",
131
"d =\n",
132
"\n",
133
" 1.00000\n",
134
" 0.50000\n",
135
"\n",
136
"y =\n",
137
"\n",
138
" 1\n",
139
" 1\n",
140
"\n"
141
]
142
}
143
],
144
"source": [
145
"A=[2,1;1,3]\n",
146
"L=[1,0;0.5,1]\n",
147
"U=[2,1;0,2.5]\n",
148
"L*U\n",
149
"\n",
150
"d=[1;0.5]\n",
151
"y=L*d"
152
]
153
},
154
{
155
"cell_type": "markdown",
156
157
"source": [
158
"## Pivoting for LU factorization\n",
159
"\n",
160
"LU factorization uses the same method as Gauss elimination so it is also necessary to perform partial pivoting when creating the lower and upper triangular matrices. \n",
161
"\n",
162
"Matlab and Octave use pivoting in the command \n",
163
"\n",
164
"`[L,U,P]=lu(A)`\n",
165
"\n"
166
]
167
},
168
{
169
"cell_type": "code",
170
"execution_count": 4,
171
172
"collapsed": false
173
},
174
"outputs": [
175
{
176
"name": "stdout",
177
"output_type": "stream",
178
"text": [
179
"'lu' is a built-in function from the file libinterp/corefcn/lu.cc\n",
180
"\n",
181
" -- Built-in Function: [L, U] = lu (A)\n",
182
" -- Built-in Function: [L, U, P] = lu (A)\n",
183
" -- Built-in Function: [L, U, P, Q] = lu (S)\n",
184
" -- Built-in Function: [L, U, P, Q, R] = lu (S)\n",
185
" -- Built-in Function: [...] = lu (S, THRES)\n",
186
" -- Built-in Function: Y = lu (...)\n",
187
" -- Built-in Function: [...] = lu (..., \"vector\")\n",
188
" Compute the LU decomposition of A.\n",
189
"\n",
190
" If A is full subroutines from LAPACK are used and if A is sparse\n",
191
" then UMFPACK is used.\n",
192
"\n",
193
" The result is returned in a permuted form, according to the\n",
194
" optional return value P. For example, given the matrix 'a = [1, 2;\n",
195
" 3, 4]',\n",
196
"\n",
197
" [l, u, p] = lu (A)\n",
198
"\n",
199
" returns\n",
200
"\n",
201
" l =\n",
202
"\n",
203
" 1.00000 0.00000\n",
204
" 0.33333 1.00000\n",
205
"\n",
206
" u =\n",
207
"\n",
208
" 3.00000 4.00000\n",
209
" 0.00000 0.66667\n",
210
"\n",
211
" p =\n",
212
"\n",
213
" 0 1\n",
214
" 1 0\n",
215
"\n",
216
" The matrix is not required to be square.\n",
217
"\n",
218
" When called with two or three output arguments and a spare input\n",
219
" matrix, 'lu' does not attempt to perform sparsity preserving column\n",
220
" permutations. Called with a fourth output argument, the sparsity\n",
221
" preserving column transformation Q is returned, such that 'P * A *\n",
222
" Q = L * U'.\n",
223
"\n",
224
" Called with a fifth output argument and a sparse input matrix, 'lu'\n",
225
" attempts to use a scaling factor R on the input matrix such that 'P\n",
226
" * (R \\ A) * Q = L * U'. This typically leads to a sparser and more\n",
227
" stable factorization.\n",
228
"\n",
229
" An additional input argument THRES, that defines the pivoting\n",
230
" threshold can be given. THRES can be a scalar, in which case it\n",
231
" defines the UMFPACK pivoting tolerance for both symmetric and\n",
232
" unsymmetric cases. If THRES is a 2-element vector, then the first\n",
233
" element defines the pivoting tolerance for the unsymmetric UMFPACK\n",
234
" pivoting strategy and the second for the symmetric strategy. By\n",
235
" default, the values defined by 'spparms' are used ([0.1, 0.001]).\n",
236
"\n",
237
" Given the string argument \"vector\", 'lu' returns the values of P\n",
238
" and Q as vector values, such that for full matrix, 'A (P,:) = L *\n",
239
" U', and 'R(P,:) * A (:, Q) = L * U'.\n",
240
"\n",
241
" With two output arguments, returns the permuted forms of the upper\n",
242
" and lower triangular matrices, such that 'A = L * U'. With one\n",
243
" output argument Y, then the matrix returned by the LAPACK routines\n",
244
" is returned. If the input matrix is sparse then the matrix L is\n",
245
" embedded into U to give a return value similar to the full case.\n",
246
" For both full and sparse matrices, 'lu' loses the permutation\n",
247
" information.\n",
248
"\n",
249
250
"\n",
251
"Additional help for built-in functions and operators is\n",
252
"available in the online version of the manual. Use the command\n",
253
"'doc <topic>' to search the manual index.\n",
254
"\n",
255
"Help and information about Octave is also available on the WWW\n",
256
"at http://www.octave.org and via the help@octave.org\n",
257
"mailing list.\n"
258
]
259
}
260
],
261
"source": [
262
"help lu"
263
]
264
},
265
{
266
"cell_type": "code",
267
"execution_count": 22,
268
269
"collapsed": false
270
},
271
"outputs": [
272
{
273
"data": {
274
"image/svg+xml": [
275
276
"\n",
277
"<title>Gnuplot</title>\n",
278
"<desc>Produced by GNUPLOT 5.0 patchlevel 3 </desc>\n",
279
"\n",
280
"<g id=\"gnuplot_canvas\">\n",
281
"\n",
282
"<rect fill=\"none\" height=\"420\" width=\"560\" x=\"0\" y=\"0\"/>\n",
283
"<defs>\n",
284
"\n",
285
"\t<circle id=\"gpDot\" r=\"0.5\" stroke-width=\"0.5\"/>\n",
286
"\t<path d=\"M-1,0 h2 M0,-1 v2\" id=\"gpPt0\" stroke=\"currentColor\" stroke-width=\"0.222\"/>\n",
287
"\t<path d=\"M-1,-1 L1,1 M1,-1 L-1,1\" id=\"gpPt1\" stroke=\"currentColor\" stroke-width=\"0.222\"/>\n",
288
"\t<path d=\"M-1,0 L1,0 M0,-1 L0,1 M-1,-1 L1,1 M-1,1 L1,-1\" id=\"gpPt2\" stroke=\"currentColor\" stroke-width=\"0.222\"/>\n",
289
"\t<rect height=\"2\" id=\"gpPt3\" stroke=\"currentColor\" stroke-width=\"0.222\" width=\"2\" x=\"-1\" y=\"-1\"/>\n",
290
"\t<rect fill=\"currentColor\" height=\"2\" id=\"gpPt4\" stroke=\"currentColor\" stroke-width=\"0.222\" width=\"2\" x=\"-1\" y=\"-1\"/>\n",
291
"\t<circle cx=\"0\" cy=\"0\" id=\"gpPt5\" r=\"1\" stroke=\"currentColor\" stroke-width=\"0.222\"/>\n",
292
293
"\t<path d=\"M0,-1.33 L-1.33,0.67 L1.33,0.67 z\" id=\"gpPt7\" stroke=\"currentColor\" stroke-width=\"0.222\"/>\n",
294
295
296
297
298
299
"\t<path d=\"M0,1.330 L1.265,0.411 L0.782,-1.067 L-0.782,-1.076 L-1.265,0.411 z\" id=\"gpPt13\" stroke=\"currentColor\" stroke-width=\"0.222\"/>\n",
300
301
"\t<filter filterUnits=\"objectBoundingBox\" height=\"1\" id=\"textbox\" width=\"1\" x=\"0\" y=\"0\">\n",
302
"\t <feFlood flood-color=\"white\" flood-opacity=\"1\" result=\"bgnd\"/>\n",
303
"\t <feComposite in=\"SourceGraphic\" in2=\"bgnd\" operator=\"atop\"/>\n",
304
"\t</filter>\n",
305
"\t<filter filterUnits=\"objectBoundingBox\" height=\"1\" id=\"greybox\" width=\"1\" x=\"0\" y=\"0\">\n",
306
"\t <feFlood flood-color=\"lightgrey\" flood-opacity=\"1\" result=\"grey\"/>\n",
307
"\t <feComposite in=\"SourceGraphic\" in2=\"grey\" operator=\"atop\"/>\n",
308
"\t</filter>\n",
309
"</defs>\n",
310
"<g color=\"white\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"1.00\">\n",
311
"</g>\n",
312
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"1.00\">\n",
313
"\t<g shape-rendering=\"crispEdges\" stroke=\"none\">\n",
314
"\t\t<polygon fill=\"rgb(255, 255, 255)\" points=\"78.8,384.0 534.9,384.0 534.9,16.8 78.8,16.8 \"/>\n",
315
"\t</g>\n",
316
"</g>\n",
317
"<g color=\"black\" fill=\"none\" stroke=\"rgb(255, 255, 255)\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
318
"</g>\n",
319
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
320
"\t<path d=\"M78.8,384.0 L91.3,384.0 M535.0,384.0 L522.5,384.0 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"end\" transform=\"translate(70.5,390.0)\">\n",
321
"\t\t<text><tspan font-family=\"{}\">0</tspan></text>\n",
322
"\t</g>\n",
323
"</g>\n",
324
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
325
"\t<path d=\"M78.8,322.8 L91.3,322.8 M535.0,322.8 L522.5,322.8 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"end\" transform=\"translate(70.5,328.8)\">\n",
326
"\t\t<text><tspan font-family=\"{}\">5e-05</tspan></text>\n",
327
"\t</g>\n",
328
"</g>\n",
329
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
330
"\t<path d=\"M78.8,261.6 L91.3,261.6 M535.0,261.6 L522.5,261.6 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"end\" transform=\"translate(70.5,267.6)\">\n",
331
"\t\t<text><tspan font-family=\"{}\">0.0001</tspan></text>\n",
332
"\t</g>\n",
333
"</g>\n",
334
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
335
"\t<path d=\"M78.8,200.3 L91.3,200.3 M535.0,200.3 L522.5,200.3 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"end\" transform=\"translate(70.5,206.3)\">\n",
336
"\t\t<text><tspan font-family=\"{}\">0.00015</tspan></text>\n",
337
"\t</g>\n",
338
"</g>\n",
339
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
340
"\t<path d=\"M78.8,139.1 L91.3,139.1 M535.0,139.1 L522.5,139.1 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"end\" transform=\"translate(70.5,145.1)\">\n",
341
"\t\t<text><tspan font-family=\"{}\">0.0002</tspan></text>\n",
342
"\t</g>\n",
343
"</g>\n",
344
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
345
"\t<path d=\"M78.8,77.9 L91.3,77.9 M535.0,77.9 L522.5,77.9 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"end\" transform=\"translate(70.5,83.9)\">\n",
346
"\t\t<text><tspan font-family=\"{}\">0.00025</tspan></text>\n",
347
"\t</g>\n",
348
"</g>\n",
349
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
350
"\t<path d=\"M78.8,16.7 L91.3,16.7 M535.0,16.7 L522.5,16.7 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"end\" transform=\"translate(70.5,22.7)\">\n",
351
"\t\t<text><tspan font-family=\"{}\">0.0003</tspan></text>\n",
352
"\t</g>\n",
353
"</g>\n",
354
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
355
"\t<path d=\"M78.8,384.0 L78.8,371.5 M78.8,16.7 L78.8,29.2 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"middle\" transform=\"translate(78.8,408.0)\">\n",
356
"\t\t<text><tspan font-family=\"{}\">0</tspan></text>\n",
357
"\t</g>\n",
358
"</g>\n",
359
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
360
"\t<path d=\"M170.0,384.0 L170.0,371.5 M170.0,16.7 L170.0,29.2 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"middle\" transform=\"translate(170.0,408.0)\">\n",
361
"\t\t<text><tspan font-family=\"{}\">20</tspan></text>\n",
362
"\t</g>\n",
363
"</g>\n",
364
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
365
"\t<path d=\"M261.3,384.0 L261.3,371.5 M261.3,16.7 L261.3,29.2 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"middle\" transform=\"translate(261.3,408.0)\">\n",
366
"\t\t<text><tspan font-family=\"{}\">40</tspan></text>\n",
367
"\t</g>\n",
368
"</g>\n",
369
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
370
"\t<path d=\"M352.5,384.0 L352.5,371.5 M352.5,16.7 L352.5,29.2 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"middle\" transform=\"translate(352.5,408.0)\">\n",
371
"\t\t<text><tspan font-family=\"{}\">60</tspan></text>\n",
372
"\t</g>\n",
373
"</g>\n",
374
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
375
"\t<path d=\"M443.8,384.0 L443.8,371.5 M443.8,16.7 L443.8,29.2 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"middle\" transform=\"translate(443.8,408.0)\">\n",
376
"\t\t<text><tspan font-family=\"{}\">80</tspan></text>\n",
377
"\t</g>\n",
378
"</g>\n",
379
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
380
"\t<path d=\"M535.0,384.0 L535.0,371.5 M535.0,16.7 L535.0,29.2 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"middle\" transform=\"translate(535.0,408.0)\">\n",
381
"\t\t<text><tspan font-family=\"{}\">100</tspan></text>\n",
382
"\t</g>\n",
383
"</g>\n",
384
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
385
"</g>\n",
386
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
387
"\t<path d=\"M78.8,16.7 L78.8,384.0 L535.0,384.0 L535.0,16.7 L78.8,16.7 Z \" stroke=\"black\"/></g>\n",
388
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
389
"</g>\n",
390
"<g color=\"black\" fill=\"none\" stroke=\"black\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"1.00\">\n",
391
"</g>\n",
392
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"1.00\">\n",
393
"\t<path d=\"M349.7,121.7 L349.7,25.7 L526.7,25.7 L526.7,121.7 L349.7,121.7 Z \" stroke=\"black\"/></g>\n",
394
"\t<g id=\"gnuplot_plot_1a\"><title>LU decomp</title>\n",
395
"<g color=\"white\" fill=\"none\" stroke=\"black\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"3.00\">\n",
396
"</g>\n",
397
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"3.00\">\n",
398
"\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"end\" transform=\"translate(526.7,55.7)\">\n",
399
"\t\t<text><tspan font-family=\"{}\">LU decomp</tspan></text>\n",
400
"\t</g>\n",
401
"</g>\n",
402
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"3.00\">\n",
403
"\t<path d=\"M360.9,49.7 L414.7,49.7 M83.4,330.3 L87.9,244.5 L92.5,289.7 L97.0,292.1 L101.6,288.5 L106.2,289.7 L110.7,287.1 L115.3,283.6 L119.9,279.8 L124.4,276.3 L129.0,276.3 L133.5,265.5 L138.1,269.0 L142.7,266.4 L147.2,261.7 L151.8,261.7 L156.4,271.3 L160.9,269.9 L165.5,274.8 L170.0,272.5 L174.6,268.7 L179.2,212.7 L183.7,256.7 L188.3,261.4 L192.9,260.2 L197.4,234.5 L202.0,270.2 L206.5,269.0 L211.1,249.4 L215.7,265.2 L220.2,303.1 L224.8,298.2 L229.3,300.8 L233.9,299.6 L238.5,292.3 L243.0,299.3 L247.6,298.2 L252.2,294.7 L256.7,284.8 L261.3,282.4 L265.8,306.9 L270.4,297.0 L275.0,303.1 L279.5,306.9 L284.1,299.3 L288.7,283.6 L293.2,294.4 L297.8,295.8 L302.3,292.3 L306.9,292.1 L311.5,283.6 L316.0,287.4 L320.6,286.2 L325.1,272.8 L329.7,293.2 L334.3,292.1 L338.8,294.7 L343.4,293.2 L348.0,292.1 L352.5,302.0 L357.1,302.0 L361.6,302.0 L366.2,299.6 L370.8,293.2 L375.3,277.5 L379.9,294.7 L384.5,292.1 L389.0,289.7 L393.6,289.7 L398.1,288.5 L402.7,286.2 L407.3,288.5 L411.8,234.5 L416.4,285.9 L421.0,282.4 L425.5,241.8 L430.1,241.8 L434.6,279.8 L439.2,266.4 L443.8,280.1 L448.3,259.1 L452.9,229.6 L457.4,277.7 L462.0,283.6 L466.6,276.3 L471.1,278.6 L475.7,277.5 L480.3,277.5 L484.8,274.8 L489.4,265.2 L493.9,270.2 L498.5,272.5 L503.1,259.1 L507.6,264.0 L512.2,264.0 L516.8,266.4 L521.3,261.4 L525.9,264.0 L530.4,247.1 L535.0,259.1 \" stroke=\"rgb( 0, 0, 255)\"/></g>\n",
404
"\t</g>\n",
405
"\t<g id=\"gnuplot_plot_2a\"><title>Octave \\\\</title>\n",
406
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"3.00\">\n",
407
"\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"end\" transform=\"translate(526.7,103.7)\">\n",
408
"\t\t<text><tspan font-family=\"{}\">Octave \\</tspan></text>\n",
409
"\t</g>\n",
410
"</g>\n",
411
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"3.00\">\n",
412
"\t<path d=\"M360.9,97.7 L414.7,97.7 M83.4,359.5 L87.9,287.4 L92.5,306.6 L97.0,311.9 L101.6,305.8 L106.2,306.9 L110.7,305.5 L115.3,303.1 L119.9,293.2 L124.4,292.3 L129.0,262.6 L133.5,278.6 L138.1,272.8 L142.7,232.2 L147.2,265.2 L151.8,260.2 L156.4,275.1 L160.9,238.3 L165.5,231.9 L170.0,240.7 L174.6,221.1 L179.2,229.9 L183.7,221.1 L188.3,218.8 L192.9,198.9 L197.4,185.8 L202.0,218.8 L206.5,221.1 L211.1,216.2 L215.7,213.8 L220.2,265.2 L224.8,265.2 L229.3,269.9 L233.9,254.1 L238.5,256.7 L243.0,256.7 L247.6,249.1 L252.2,247.1 L256.7,241.8 L261.3,204.2 L265.8,265.2 L270.4,262.9 L275.0,243.3 L279.5,264.0 L284.1,252.9 L288.7,249.4 L293.2,247.1 L297.8,245.6 L302.3,241.0 L306.9,229.9 L311.5,205.4 L316.0,229.9 L320.6,205.1 L325.1,222.3 L329.7,241.0 L334.3,244.2 L338.8,238.3 L343.4,227.2 L348.0,218.5 L352.5,250.6 L357.1,250.6 L361.6,246.8 L366.2,232.2 L370.8,240.7 L375.3,238.3 L379.9,220.0 L384.5,207.7 L389.0,222.6 L393.6,217.6 L398.1,215.0 L402.7,212.7 L407.3,213.8 L411.8,223.5 L416.4,208.9 L421.0,205.4 L425.5,196.6 L430.1,184.6 L434.6,190.5 L439.2,195.4 L443.8,185.8 L448.3,177.0 L452.9,183.2 L457.4,180.8 L462.0,188.1 L466.6,164.8 L471.1,164.8 L475.7,175.9 L480.3,175.9 L484.8,169.7 L489.4,165.9 L493.9,161.0 L498.5,75.5 L503.1,150.2 L507.6,135.3 L512.2,150.2 L516.8,145.5 L521.3,147.6 L525.9,134.1 L530.4,126.8 L535.0,135.3 \" stroke=\"rgb( 0, 128, 0)\"/></g>\n",
413
"\t</g>\n",
414
"<g color=\"white\" fill=\"none\" stroke=\"rgb( 0, 128, 0)\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"2.00\">\n",
415
"</g>\n",
416
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"2.00\">\n",
417
"</g>\n",
418
"<g color=\"black\" fill=\"none\" stroke=\"black\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
419
"</g>\n",
420
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
421
"</g>\n",
422
"</g>\n",
423
"</svg>"
424
],
425
"text/plain": [
426
"<IPython.core.display.SVG object>"
427
]
428
},
429
430
"output_type": "display_data"
431
}
432
],
433
"source": [
434
"% time LU solution vs backslash\n",
435
"t_lu=zeros(100,1);\n",
436
"t_bs=zeros(100,1);\n",
437
"for N=1:100\n",
438
" A=rand(N,N);\n",
439
" y=rand(N,1);\n",
440
" [L,U,P]=lu(A);\n",
441
"\n",
442
" tic; d=L\\y; x=U\\d; t_lu(N)=toc;\n",
443
"\n",
444
" tic; x=A\\y; t_bs(N)=toc;\n",
445
"end\n",
446
"plot([1:100],t_lu,[1:100],t_bs) \n",
447
"legend('LU decomp','Octave \\\\')"
448
]
449
},
450
{
451
"cell_type": "markdown",
452
453
"collapsed": true
454
},
455
"source": [
456
"Consider the problem again from the intro to Linear Algebra, 4 masses are connected in series to 4 springs with K=10 N/m. What are the final positions of the masses? \n",
457
"\n",
458
"![Springs-masses](../lecture_09/mass_springs.svg)\n",
459
"\n",
460
"The masses haves the following amounts, 1, 2, 3, and 4 kg for masses 1-4. Using a FBD for each mass:\n",
461
"\n",
462
"\$m_{1}g+k(x_{2}-x_{1})-kx_{1}=0\$\n",
463
"\n",
464
"\$m_{2}g+k(x_{3}-x_{2})-k(x_{2}-x_{1})=0\$\n",
465
"\n",
466
"\$m_{3}g+k(x_{4}-x_{3})-k(x_{3}-x_{2})=0\$\n",
467
"\n",
468
"\$m_{4}g-k(x_{4}-x_{3})=0\$\n",
469
"\n",
470
"in matrix form:\n",
471
"\n",
472
"\$\\left[ \\begin{array}{cccc}\n",
473
"2k & -k & 0 & 0 \\\\\n",
474
"-k & 2k & -k & 0 \\\\\n",
475
"0 & -k & 2k & -k \\\\\n",
476
"0 & 0 & -k & k \\end{array} \\right]\n",
477
"\\left[ \\begin{array}{c}\n",
478
"x_{1} \\\\\n",
479
"x_{2} \\\\\n",
480
"x_{3} \\\\\n",
481
"x_{4} \\end{array} \\right]=\n",
482
"\\left[ \\begin{array}{c}\n",
483
"m_{1}g \\\\\n",
484
"m_{2}g \\\\\n",
485
"m_{3}g \\\\\n",
486
"m_{4}g \\end{array} \\right]\$"
487
]
488
},
489
{
490
"cell_type": "code",
491
"execution_count": 24,
492
493
"collapsed": false
494
},
495
"outputs": [
496
{
497
"name": "stdout",
498
"output_type": "stream",
499
"text": [
500
"K =\n",
501
"\n",
502
" 20 -10 0 0\n",
503
" -10 20 -10 0\n",
504
" 0 -10 20 -10\n",
505
" 0 0 -10 10\n",
506
"\n",
507
"y =\n",
508
"\n",
509
" 9.8100\n",
510
" 19.6200\n",
511
" 29.4300\n",
512
" 39.2400\n",
513
"\n"
514
]
515
}
516
],
517
"source": [
518
"k=10; % N/m\n",
519
"m1=1; % kg\n",
520
"m2=2;\n",
521
"m3=3;\n",
522
"m4=4;\n",
523
"g=9.81; % m/s^2\n",
524
"K=[2*k -k 0 0; -k 2*k -k 0; 0 -k 2*k -k; 0 0 -k k]\n",
525
"y=[m1*g;m2*g;m3*g;m4*g]"
526
]
527
},
528
{
529
"cell_type": "markdown",
530
531
"source": [
532
"This matrix, K, is symmetric. \n",
533
"\n",
534
"`K(i,j)==K(j,i)`\n",
535
"\n",
536
"Now we can use,\n",
537
"\n",
538
"## Cholesky Factorization\n",
539
"\n",
540
"We can decompose the matrix, K into two matrices, \$U\$ and \$U^{T}\$, where\n",
541
"\n",
542
"\$K=U^{T}U\$\n",
543
"\n",
544
"each of the components of U can be calculated with the following equations:\n",
545
"\n",
546
"\$u_{ii}=\\sqrt{a_{ii}-\\sum_{k=1}^{i-1}u_{ki}^{2}}\$\n",
547
"\n",
548
"\$u_{ij}=\\frac{a_{ij}-\\sum_{k=1}^{i-1}u_{ki}u_{kj}}{u_{ii}}\$\n",
549
"\n",
550
"so for K"
551
]
552
},
553
{
554
"cell_type": "code",
555
"execution_count": 25,
556
557
"collapsed": false
558
},
559
"outputs": [
560
{
561
"name": "stdout",
562
"output_type": "stream",
563
"text": [
564
"K =\n",
565
"\n",
566
" 20 -10 0 0\n",
567
" -10 20 -10 0\n",
568
" 0 -10 20 -10\n",
569
" 0 0 -10 10\n",
570
"\n"
571
]
572
}
573
],
574
"source": [
575
"K"
576
]
577
},
578
{
579
"cell_type": "code",
580
"execution_count": 26,
581
582
"collapsed": false
583
},
584
"outputs": [
585
{
586
"name": "stdout",
587
"output_type": "stream",
588
"text": [
589
"u11 = 4.4721\n",
590
"u12 = -2.2361\n",
591
"u13 = 0\n",
592
"u14 = 0\n",
593
"u22 = 3.8730\n",
594
"u23 = -2.5820\n",
595
"u24 = 0\n",
596
"u33 = 3.6515\n",
597
"u34 = -2.7386\n",
598
"u44 = 1.5811\n",
599
"U =\n",
600
"\n",
601
" 4.47214 -2.23607 0.00000 0.00000\n",
602
" 0.00000 3.87298 -2.58199 0.00000\n",
603
" 0.00000 0.00000 3.65148 -2.73861\n",
604
" 0.00000 0.00000 0.00000 1.58114\n",
605
"\n"
606
]
607
}
608
],
609
"source": [
610
"u11=sqrt(K(1,1))\n",
611
"u12=(K(1,2))/u11\n",
612
"u13=(K(1,3))/u11\n",
613
"u14=(K(1,4))/u11\n",
614
"u22=sqrt(K(2,2)-u12^2)\n",
615
"u23=(K(2,3)-u12*u13)/u22\n",
616
"u24=(K(2,4)-u12*u14)/u22\n",
617
"u33=sqrt(K(3,3)-u13^2-u23^2)\n",
618
"u34=(K(3,4)-u13*u14-u23*u24)/u33\n",
619
"u44=sqrt(K(4,4)-u14^2-u24^2-u34^2)\n",
620
"U=[u11,u12,u13,u14;0,u22,u23,u24;0,0,u33,u34;0,0,0,u44]"
621
]
622
},
623
{
624
"cell_type": "code",
625
"execution_count": 27,
626
627
"collapsed": false
628
},
629
"outputs": [
630
{
631
"name": "stdout",
632
"output_type": "stream",
633
"text": [
634
"ans =\n",
635
"\n",
636
" 20.00000 -10.00000 0.00000 0.00000\n",
637
" -10.00000 20.00000 -10.00000 0.00000\n",
638
" 0.00000 -10.00000 20.00000 -10.00000\n",
639
" 0.00000 0.00000 -10.00000 10.00000\n",
640
"\n"
641
]
642
}
643
],
644
"source": [
645
"U'*U"
646
]
647
},
648
{
649
"cell_type": "code",
650
"execution_count": 37,
651
652
"collapsed": false
653
},
654
"outputs": [
655
{
656
"name": "stdout",
657
"output_type": "stream",
658
"text": [
659
"average time spent for Cholesky factored solution = 1.623154e-05+/-1.166726e-05\n",
660
"average time spent for backslash solution = 1.675844e-05+/-1.187234e-05\n"
661
]
662
}
663
],
664
"source": [
665
"% time solution for Cholesky vs backslash\n",
666
"t_chol=zeros(1000,1);\n",
667
"t_bs=zeros(1000,1);\n",
668
"for i=1:1000\n",
669
" tic; d=U'*y; x=U\\d; t_chol(i)=toc;\n",
670
" tic; x=K\\y; t_bs(i)=toc;\n",
671
"end\n",
672
"fprintf('average time spent for Cholesky factored solution = %e+/-%e',mean(t_chol),std(t_chol))\n",
673
"\n",
674
"fprintf('average time spent for backslash solution = %e+/-%e',mean(t_bs),std(t_bs))"
675
]
676
}
677
],
678
679
"kernelspec": {
680
"display_name": "Octave",
681
"language": "octave",
682
"name": "octave"
683
},
684
"language_info": {
685
"file_extension": ".m",
686
687
{
688
"text": "MetaKernel Magics",
689
690
}
691
],
692
"mimetype": "text/x-octave",
693
"name": "octave",
694
"version": "0.19.14"
695
}
696
},
697
"nbformat": 4,
698
"nbformat_minor": 2
699
}
You can’t perform that action at this time.