-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathindex.html
46 lines (46 loc) · 20.6 KB
/
index.html
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
<!DOCTYPE html><html lang="en"><head><meta http-equiv="content-type" content="text/html; charset=utf-8"><meta content="width=device-width, initial-scale=1.0, maximum-scale=1.0, user-scalable=0" name="viewport"><meta content="yes" name="apple-mobile-web-app-capable"><meta content="black-translucent" name="apple-mobile-web-app-status-bar-style"><meta content="telephone=no" name="format-detection"><meta name="description"><title>迭代法解线性方程组 | Welcome To Oa</title><link rel="stylesheet" type="text/css" href="/css/style.css?v=0.0.0"><link rel="stylesheet" type="text/css" href="/css/donate.css"><link rel="stylesheet" type="text/css" href="//cdn.bootcss.com/normalize/6.0.0/normalize.min.css"><link rel="stylesheet" type="text/css" href="//cdn.bootcss.com/pure/0.6.2/pure-min.css"><link rel="stylesheet" type="text/css" href="//cdn.bootcss.com/pure/0.6.2/grids-responsive-min.css"><link rel="stylesheet" href="//cdn.bootcss.com/font-awesome/4.7.0/css/font-awesome.min.css"><script type="text/javascript" src="//cdn.bootcss.com/jquery/3.2.1/jquery.min.js"></script><link rel="Shortcut Icon" type="image/x-icon" href="/favicon.ico"><link rel="apple-touch-icon" href="/apple-touch-icon.png"><link rel="apple-touch-icon-precomposed" href="/apple-touch-icon.png"></head><body><div class="body_container"><div id="header"><div class="site-name"><h1 class="hidden">迭代法解线性方程组</h1><a id="logo" href="/.">Welcome To Oa</a><p class="description"></p></div><div id="nav-menu"><a href="/." class="current"><i class="fa fa-home"> Home</i></a><a href="/archives/"><i class="fa fa-archive"> Archive</i></a><a href="/about/"><i class="fa fa-user"> About</i></a><a href="/guestbook/"><i class="fa fa-comments"> Guestbook</i></a><a href="/atom.xml"><i class="fa fa-rss"> RSS</i></a></div></div><div id="layout" class="pure-g"><div class="pure-u-1 pure-u-md-3-4"><div class="content_container"><div class="post"><h1 class="post-title">迭代法解线性方程组</h1><div class="post-meta">Oct 28, 2017<script src="https://dn-lbstatics.qbox.me/busuanzi/2.3/busuanzi.pure.mini.js" async></script><span id="busuanzi_container_page_pv"> | <span id="busuanzi_value_page_pv"></span><span> Hits</span></span></div><div class="clear"><div id="toc" class="toc-article"><div class="toc-title">Contents</div><ol class="toc"><li class="toc-item toc-level-1"><a class="toc-link" href="#undefined"><span class="toc-number">1.</span> <span class="toc-text">要求</span></a></li><li class="toc-item toc-level-1"><a class="toc-link" href="#undefined"><span class="toc-number">2.</span> <span class="toc-text">分析</span></a><ol class="toc-child"><li class="toc-item toc-level-2"><a class="toc-link" href="#undefined"><span class="toc-number">2.1.</span> <span class="toc-text">Jacobi迭代法</span></a></li><li class="toc-item toc-level-2"><a class="toc-link" href="#undefined"><span class="toc-number">2.2.</span> <span class="toc-text">Gauss-Seidel迭代法</span></a></li><li class="toc-item toc-level-2"><a class="toc-link" href="#undefined"><span class="toc-number">2.3.</span> <span class="toc-text">SOR</span></a></li></ol></li><li class="toc-item toc-level-1"><a class="toc-link" href="#undefined"><span class="toc-number">3.</span> <span class="toc-text">实现</span></a><ol class="toc-child"><li class="toc-item toc-level-2"><a class="toc-link" href="#undefined"><span class="toc-number">3.1.</span> <span class="toc-text">Basic</span></a></li><li class="toc-item toc-level-2"><a class="toc-link" href="#undefined"><span class="toc-number">3.2.</span> <span class="toc-text">Jacobi</span></a></li><li class="toc-item toc-level-2"><a class="toc-link" href="#undefined"><span class="toc-number">3.3.</span> <span class="toc-text">Gauss-Seidel</span></a></li><li class="toc-item toc-level-2"><a class="toc-link" href="#undefined"><span class="toc-number">3.4.</span> <span class="toc-text">SOR</span></a></li><li class="toc-item toc-level-2"><a class="toc-link" href="#undefined"><span class="toc-number">3.5.</span> <span class="toc-text">Main</span></a></li></ol></li><li class="toc-item toc-level-1"><a class="toc-link" href="#undefined"><span class="toc-number">4.</span> <span class="toc-text">结果分析</span></a></li></ol></div></div><div class="post-content"><h1>要求</h1>
<p>分别用①Jacobi，②Gauss-Seidel，③SOR解线性方程组</p>
<h1>分析</h1>
<h2>Jacobi迭代法</h2>
<p>输入：Ax=b
所以$$a_{ii}x_i^{(k+1)}=b_i-\sum_{j=1&&j!=i}^na_{ij}x_j^{(k)} $$
将系数矩阵A表示为$$A=D-L-U$$
其中D为主对角线元素为A，其他为0的矩阵，-L，-U分别为下三角等于A，上三角等于A，其他元素为0的矩阵。又因为D可逆\
因此我们可得到$$X^{(k+1)}=D^{-1}(L+U)X^{(k)}+D^{-1}b$$</p>
<h2>Gauss-Seidel迭代法</h2>
<p>输入：Ax=b
与Jacobi类似，但我们用$x_i^{(k+1)}$代替$x_i^{(k)}$进行估计所以$$x_i^{(k+1)}=1/a_{ii}(b_i-\sum_{j=1}^{i-1}a_{ij}x_j^{(k+1)}-\sum_{j=i+1}^ia_{ij}x_j^{(k)})$$
即$$X^{(k+1)}=(D-L)^{-1}UX^{(k)}+(D-L)^{-1}b$$</p>
<h2>SOR</h2>
<p>通过添加合适的松弛因子w，使收敛速度变快
输入：Ax=b
与Gauss-Seidel类似，再通过添加合适的松弛因子w，使收敛速度变快。所以$$x_i^{(k+1)}=1/a_{ii}((1-w)a_{ii}x_i^{(k)}+w(b_i-\sum_{j=1}^{i-1}a_{ij}x_j^{(k+1)}-\sum_{j=i+1}^na_{ij}x_j^{(k)}))$$即$$X^{(k+1)}=(D-wL)^{-1}[(1-w)D+wU]X^{(k)}+w(D-wL)^{-1}b$$</p>
<h1>实现</h1>
<h2>Basic</h2>
<p><figure class="highlight python"><table><tr><td class="gutter"><pre><span class="line">1</span><br><span class="line">2</span><br><span class="line">3</span><br><span class="line">4</span><br><span class="line">5</span><br><span class="line">6</span><br><span class="line">7</span><br><span class="line">8</span><br><span class="line">9</span><br><span class="line">10</span><br><span class="line">11</span><br><span class="line">12</span><br><span class="line">13</span><br><span class="line">14</span><br></pre></td><td class="code"><pre><span class="line"><span class="comment">#decomposite A to D,L,U</span></span><br><span class="line"><span class="function"><span class="keyword">def</span> <span class="title">getDLU</span><span class="params">(A)</span>:</span></span><br><span class="line"> D=np.zeros(A.shape)</span><br><span class="line"> L=np.zeros(A.shape)</span><br><span class="line"> U=np.zeros(A.shape)</span><br><span class="line"> <span class="keyword">for</span> i <span class="keyword">in</span> range(<span class="number">0</span>,A.shape[<span class="number">0</span>]):</span><br><span class="line"> <span class="keyword">for</span> j <span class="keyword">in</span> range(<span class="number">0</span>,A.shape[<span class="number">1</span>]):</span><br><span class="line"> <span class="keyword">if</span>(i==j):</span><br><span class="line"> D[i,j]=A[i,j]</span><br><span class="line"> <span class="keyword">elif</span>(i>j):</span><br><span class="line"> L[i,j]=-A[i,j]</span><br><span class="line"> <span class="keyword">else</span>:</span><br><span class="line"> U[i,j]=-A[i,j]</span><br><span class="line"> <span class="keyword">return</span> D,L,U</span><br></pre></td></tr></table></figure></p>
<h2>Jacobi</h2>
<p><figure class="highlight python"><table><tr><td class="gutter"><pre><span class="line">1</span><br><span class="line">2</span><br><span class="line">3</span><br><span class="line">4</span><br><span class="line">5</span><br><span class="line">6</span><br><span class="line">7</span><br><span class="line">8</span><br><span class="line">9</span><br><span class="line">10</span><br><span class="line">11</span><br><span class="line">12</span><br><span class="line">13</span><br></pre></td><td class="code"><pre><span class="line"><span class="comment">#use Jacobi iteration get x_result</span></span><br><span class="line"><span class="function"><span class="keyword">def</span> <span class="title">Jacobi</span><span class="params">(A,x,b)</span>:</span></span><br><span class="line"> D,L,U=getDLU(A)</span><br><span class="line"> D_inv=np.linalg.inv(D)</span><br><span class="line"> x_result=np.dot(np.dot(D_inv,(L+U)),x)+np.dot(D_inv,b)</span><br><span class="line"> <span class="comment">#calculate the rate</span></span><br><span class="line"> rate=np.sqrt(np.sum((np.dot(A,x_result)-b)**<span class="number">2</span>))/np.sqrt(np.sum(b**<span class="number">2</span>))</span><br><span class="line"> <span class="comment"># print rate</span></span><br><span class="line"> <span class="keyword">while</span> rate><span class="number">10</span>**<span class="number">-6</span>:</span><br><span class="line"> x=x_result</span><br><span class="line"> x_result=np.dot(np.dot(D_inv,(L+U)),x)+np.dot(D_inv,b)</span><br><span class="line"> rate=np.sqrt(np.sum((np.dot(A,x_result)-b)**<span class="number">2</span>))/np.sqrt(np.sum(b**<span class="number">2</span>))</span><br><span class="line"> <span class="keyword">return</span> x_result</span><br></pre></td></tr></table></figure></p>
<h2>Gauss-Seidel</h2>
<p><figure class="highlight python"><table><tr><td class="gutter"><pre><span class="line">1</span><br><span class="line">2</span><br><span class="line">3</span><br><span class="line">4</span><br><span class="line">5</span><br><span class="line">6</span><br><span class="line">7</span><br><span class="line">8</span><br><span class="line">9</span><br><span class="line">10</span><br><span class="line">11</span><br><span class="line">12</span><br><span class="line">13</span><br></pre></td><td class="code"><pre><span class="line"><span class="comment">#use GaussSeidel iteration get x_result</span></span><br><span class="line"><span class="function"><span class="keyword">def</span> <span class="title">GaussSeidel</span><span class="params">(A,x,b)</span>:</span></span><br><span class="line"> D,L,U=getDLU(A)</span><br><span class="line"> DL_inv=np.linalg.inv(D-L)</span><br><span class="line"> x_result=np.dot(np.dot(DL_inv,U),x)+np.dot(DL_inv,b)</span><br><span class="line"> <span class="comment">#calculate the rate</span></span><br><span class="line"> rate=np.sqrt(np.sum((np.dot(A,x_result)-b)**<span class="number">2</span>))/np.sqrt(np.sum(b**<span class="number">2</span>))</span><br><span class="line"> <span class="keyword">while</span> rate><span class="number">10</span>**<span class="number">-6</span>:</span><br><span class="line"> x=x_result</span><br><span class="line"> x_result=np.dot(np.dot(DL_inv,U),x)+np.dot(DL_inv,b)</span><br><span class="line"> rate=np.sqrt(np.sum((np.dot(A,x_result)-b)**<span class="number">2</span>))/np.sqrt(np.sum(b**<span class="number">2</span>))</span><br><span class="line"> <span class="comment"># print rate</span></span><br><span class="line"> <span class="keyword">return</span> x_result</span><br></pre></td></tr></table></figure></p>
<h2>SOR</h2>
<p><figure class="highlight python"><table><tr><td class="gutter"><pre><span class="line">1</span><br><span class="line">2</span><br><span class="line">3</span><br><span class="line">4</span><br><span class="line">5</span><br><span class="line">6</span><br><span class="line">7</span><br><span class="line">8</span><br><span class="line">9</span><br><span class="line">10</span><br><span class="line">11</span><br><span class="line">12</span><br></pre></td><td class="code"><pre><span class="line"><span class="comment">#use SOR iteration get x_result</span></span><br><span class="line"><span class="function"><span class="keyword">def</span> <span class="title">SOR</span><span class="params">(A,x,b,w)</span>:</span></span><br><span class="line"> D,L,U=getDLU(A)</span><br><span class="line"> DLw_inv=np.linalg.inv(D-w*L)</span><br><span class="line"> x_result=np.dot(np.dot(DLw_inv,(<span class="number">1</span>-w)*D+w*U),x)+np.dot(w*DLw_inv,b)</span><br><span class="line"> <span class="comment">#calculate the rate</span></span><br><span class="line"> rate=np.sqrt(np.sum((np.dot(A,x_result)-b)**<span class="number">2</span>))/np.sqrt(np.sum(b**<span class="number">2</span>))</span><br><span class="line"> <span class="keyword">while</span> rate><span class="number">10</span>**<span class="number">-6</span>:</span><br><span class="line"> x=x_result</span><br><span class="line"> x_result=np.dot(np.dot(DLw_inv,(<span class="number">1</span>-w)*D+w*U),x)+np.dot(w*DLw_inv,b)</span><br><span class="line"> rate=np.sqrt(np.sum((np.dot(A,x_result)-b)**<span class="number">2</span>))/np.sqrt(np.sum(b**<span class="number">2</span>))</span><br><span class="line"> <span class="keyword">return</span> x_result</span><br></pre></td></tr></table></figure></p>
<h2>Main</h2>
<p><figure class="highlight python"><table><tr><td class="gutter"><pre><span class="line">1</span><br><span class="line">2</span><br><span class="line">3</span><br><span class="line">4</span><br><span class="line">5</span><br><span class="line">6</span><br><span class="line">7</span><br><span class="line">8</span><br><span class="line">9</span><br><span class="line">10</span><br><span class="line">11</span><br><span class="line">12</span><br><span class="line">13</span><br><span class="line">14</span><br><span class="line">15</span><br><span class="line">16</span><br><span class="line">17</span><br><span class="line">18</span><br><span class="line">19</span><br><span class="line">20</span><br></pre></td><td class="code"><pre><span class="line"><span class="keyword">if</span> __name__==<span class="string">'__main__'</span>:</span><br><span class="line"> A=generateA()</span><br><span class="line"> b=generateb()</span><br><span class="line"> <span class="comment"># print np.dot(np.linalg.inv(A),b)</span></span><br><span class="line"> x=np.zeros(<span class="number">100</span>).T</span><br><span class="line"> time0=time.time()</span><br><span class="line"> x_result=Jacobi(A,x,b)</span><br><span class="line"> time1=time.time()-time0</span><br><span class="line"> <span class="keyword">print</span> <span class="string">"Jacobi"</span> ,time1</span><br><span class="line"> <span class="keyword">print</span> x_result</span><br><span class="line"> time0=time.time()</span><br><span class="line"> x_result=GaussSeidel(A,x,b)</span><br><span class="line"> time1=time.time()-time0</span><br><span class="line"> <span class="keyword">print</span> <span class="string">"GaussSeidel"</span>,time1</span><br><span class="line"> <span class="keyword">print</span> x_result</span><br><span class="line"> time0=time.time()</span><br><span class="line"> x_result=SOR(A,x,b,<span class="number">1.4</span>)</span><br><span class="line"> time1=time.time()-time0</span><br><span class="line"> <span class="keyword">print</span> <span class="string">"SOR"</span>,time1</span><br><span class="line"> <span class="keyword">print</span> x_result</span><br></pre></td></tr></table></figure></p>
<h1>结果分析</h1>
<p>Jacobi 3.02215313911
GaussSeidel 1.40041995049
SOR 0.781527996063
都收敛到了目的值附近，但3个收敛速度不同，其中SOR收敛最快，应对大型稀疏矩阵更加适用</p>
</div><script type="text/javascript" src="/js/share.js?v=0.0.0" async></script><a data-url="https://jarvis.world/迭代法解线性方程组/" data-id="cjv2m69ca001e1xihriih3iy4" class="article-share-link">Share</a><div class="tags"><a href="/tags/数值计算与优化/">数值计算与优化</a></div><div class="post-nav"><a href="/DL初级数学基础/" class="pre">DL初级数学基础整理</a><a href="/共轭梯度法-QR分解法/" class="next">共轭梯度法&QR分解法</a></div><div id="container"></div><link rel="stylesheet" href="https://imsun.github.io/gitment/style/default.css"><script src="https://imsun.github.io/gitment/dist/gitment.browser.js"></script><script>var gitment = new Gitment({
id: '迭代法解线性方程组',
owner: 'Jarvis-K',
repo: 'Jarvis-K.github.io',
oauth: {
client_id: 'f5c6bbd5b2ae20c2f136',
client_secret: 'fed7dc675f5293f90d56d4055654fb89db9003a8',
},
})
gitment.render('container')</script></div></div></div><div class="pure-u-1-4 hidden_mid_and_down"><div id="sidebar"><div class="widget"><form action="//www.google.com/search" method="get" accept-charset="utf-8" target="_blank" class="search-form"><input type="text" name="q" maxlength="20" placeholder="Search"/><input type="hidden" name="sitesearch" value="https://jarvis.world"/></form></div><div class="widget"><div class="widget-title"><i class="fa fa-folder-o"> Categories</i></div></div><div class="widget"><div class="widget-title"><i class="fa fa-star-o"> Tags</i></div><div class="tagcloud"><a href="/tags/util/" style="font-size: 15px;">util</a> <a href="/tags/Math/" style="font-size: 15px;">Math</a> <a href="/tags/ml-algo/" style="font-size: 15px;">ml-algo</a> <a href="/tags/gnn-papers/" style="font-size: 15px;">gnn-papers</a> <a href="/tags/installation/" style="font-size: 15px;">installation</a> <a href="/tags/Util/" style="font-size: 15px;">Util</a> <a href="/tags/rl-papers/" style="font-size: 15px;">rl-papers</a> <a href="/tags/mongodb/" style="font-size: 15px;">mongodb</a> <a href="/tags/DL/" style="font-size: 15px;">DL</a> <a href="/tags/python/" style="font-size: 15px;">python</a> <a href="/tags/mxnet/" style="font-size: 15px;">mxnet</a> <a href="/tags/gluon/" style="font-size: 15px;">gluon</a> <a href="/tags/pyqt/" style="font-size: 15px;">pyqt</a> <a href="/tags/cv/" style="font-size: 15px;">cv</a> <a href="/tags/数值计算与优化/" style="font-size: 15px;">数值计算与优化</a> <a href="/tags/ML/" style="font-size: 15px;">ML</a> <a href="/tags/DM/" style="font-size: 15px;">DM</a></div></div><div class="widget"><div class="widget-title"><i class="fa fa-file-o"> Recent</i></div><ul class="post-list"><li class="post-list-item"><a class="post-list-link" href="/Hypergraph-Convolution-and-Hyper-Attention/">Hypergraph Convolution and Hyper Attention</a></li><li class="post-list-item"><a class="post-list-link" href="/EM-Algorithm/">EM Algorithm</a></li><li class="post-list-item"><a class="post-list-link" href="/Summary-Mean-Field-Multi-Agent-Reinforcement-Learning/">Summary: Mean Field Multi-Agent Reinforcement Learning. ICML(2018)</a></li><li class="post-list-item"><a class="post-list-link" href="/Summary-ACCNet-Actor-Coordinator-Critic-Net-for-“Learning-to-Communicate”-with-Deep-Multi-agent-Reinforcement-Learning/">Summary: ACCNet: Actor-Coordinator-Critic Net for “Learning-to-Communicate” with Deep Multi-agent Reinforcement Learning</a></li><li class="post-list-item"><a class="post-list-link" href="/Summary-Emergence-of-Grounded-Compositional-Language-in-Multi-Agent-Populations/">Summary: Emergence of Grounded Compositional Language in Multi-Agent Populations(AAAI 2018)</a></li><li class="post-list-item"><a class="post-list-link" href="/Summary-Learning-Multiagent-Communication-with-Backpropagation/">Summary: Learning Multiagent Communication with Backpropagation(Nips 2016)</a></li><li class="post-list-item"><a class="post-list-link" href="/Summary-Learning-to-Communicate-with-Deep-Multi-Agent-Reinforcement-Learning/">Summary: Learning to Communicate with Deep Multi-Agent Reinforcement Learning(Nips 2016)</a></li><li class="post-list-item"><a class="post-list-link" href="/Summary-TarMAC-Targeted-Multi-Agent-Communication/">Summary: TarMAC:Targeted Multi-Agent Communication</a></li><li class="post-list-item"><a class="post-list-link" href="/Summary-Intrinsic-Social-Motivation-Via-Causal-Influence-In-Multi-Agent-RL/">Summary: Intrinsic Social Motivation Via Causal Influence In Multi-Agent RL</a></li><li class="post-list-item"><a class="post-list-link" href="/Summary-Actor-Attention-Critic-For-Multi-Agent-Reinforcement-Learning/">Summary: Actor-Attention-Critic For Multi-Agent Reinforcement Learning</a></li></ul></div><div class="widget"><div class="widget-title"><i class="fa fa-external-link"> Links</i></div><ul></ul><a href="https://cww97.cn/" title="cww" target="_blank">cww</a></div></div></div><div class="pure-u-1 pure-u-md-3-4"><div id="footer">Copyright © 2019 <a href="/." rel="nofollow">Welcome To Oa.</a> Powered by<a rel="nofollow" target="_blank" href="https://hexo.io"> Hexo.</a><a rel="nofollow" target="_blank" href="https://github.com/tufu9441/maupassant-hexo"> Theme</a> by<a rel="nofollow" target="_blank" href="https://github.com/pagecho"> Cho.</a></div></div></div><a id="rocket" href="#top" class="show"></a><script type="text/javascript" src="/js/totop.js?v=0.0.0" async></script><script type="text/javascript" src="//cdn.bootcss.com/fancybox/2.1.5/jquery.fancybox.pack.js" async></script><script type="text/javascript" src="/js/fancybox.js?v=0.0.0" async></script><link rel="stylesheet" type="text/css" href="/css/jquery.fancybox.css?v=0.0.0"><script type="text/x-mathjax-config">MathJax.Hub.Config({
tex2jax: {inlineMath: [['$','$'], ['\\(','\\)']]}
});
</script><script type="text/javascript" src="https://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-MML-AM_CHTML" async></script><script type="text/javascript" src="/js/codeblock-resizer.js?v=0.0.0"></script><script type="text/javascript" src="/js/smartresize.js?v=0.0.0"></script></div></body></html>