<html xmlns="http://www.w3.org/1999/xhtml">
<head>
<meta charset="utf-8">
<meta name="generator" content="pdf2htmlEX">
<meta http-equiv="X-UA-Compatible" content="IE=edge,chrome=1">
<link rel="stylesheet" href="https://static.pudn.com/base/css/base.min.css">
<link rel="stylesheet" href="https://static.pudn.com/base/css/fancy.min.css">
<link rel="stylesheet" href="https://static.pudn.com/prod/directory_preview_static/6273a38a40256a40ce34026f/raw.css">
<script src="https://static.pudn.com/base/js/compatibility.min.js"></script>
<script src="https://static.pudn.com/base/js/pdf2htmlEX.min.js"></script>
<script>
try{
pdf2htmlEX.defaultViewer = new pdf2htmlEX.Viewer({});
}catch(e){}
</script>
<title></title>
</head>
<body>
<div id="sidebar" style="display: none">
<div id="outline">
</div>
</div>
<div id="pf1" class="pf w0 h0" data-page-no="1"><div class="pc pc1 w0 h0"><img class="bi x0 y0 w1 h1" alt="" src="https://static.pudn.com/prod/directory_preview_static/6273a38a40256a40ce34026f/bg1.jpg"><div class="c x0 y1 w2 h2"><div class="t m0 x1 h3 y2 ff1 fs0 fc0 sc0 ls0 ws0">第三题：</div><div class="t m0 x2 h4 y3 ff2 fs1 fc0 sc1 ls0 ws0">1.<span class="_ _0"> </span><span class="ff1">算法设计方案</span></div><div class="t m0 x2 h5 y4 ff1 fs2 fc0 sc1 ls0 ws0">（<span class="ff3">1<span class="_ _1"></span></span>）<span class="ff3">.<span class="_ _1"></span></span>解非<span class="_ _1"></span>线性<span class="_ _1"></span>方<span class="_ _1"></span>程组<span class="_ _1"></span>：开<span class="_ _1"></span>始设<span class="_ _1"></span>计的<span class="_ _1"></span>时候<span class="_ _1"></span>使用<span class="_ _1"></span>简单<span class="_ _1"></span>迭代<span class="_ _1"></span>法，<span class="_ _1"></span>结果<span class="_ _1"></span>程序<span class="_ _1"></span>调<span class="_ _1"></span>试过<span class="_ _1"></span>冲</div><div class="t m0 x2 h5 y5 ff1 fs2 fc0 sc1 ls0 ws0">中发<span class="_ _1"></span>现简<span class="_ _1"></span>单迭<span class="_ _1"></span>代法<span class="_ _1"></span>不<span class="_ _1"></span>收敛<span class="_ _1"></span>，所<span class="_ _1"></span>以改<span class="_ _1"></span>用连<span class="_ _1"></span>续牛<span class="_ _1"></span>顿迭<span class="_ _1"></span>代法<span class="_ _1"></span>，迭<span class="_ _1"></span>代初<span class="_ _1"></span>始向<span class="_ _1"></span>量均<span class="_ _1"></span>选<span class="_ _1"></span>为：</div><div class="t m0 x2 h5 y6 ff3 fs2 fc0 sc1 ls0 ws0">x0=<span class="_ _2"> </span><span class="ff1">（<span class="_ _2"> </span></span>1,2,1,2<span class="_ _2"> </span><span class="ff1">）<span class="_ _2"> </span>；<span class="_ _2"> </span>对<span class="_ _2"> </span>原<span class="_ _2"> </span>题<span class="_ _2"> </span>中<span class="_ _2"> </span>给<span class="_ _3"> </span>出<span class="_ _3"> </span>的<span class="_ _4"> </span>，<span class="_ _5"> </span>，</span></div><div class="t m0 x3 h5 y7 ff1 fs2 fc0 sc1 ls0 ws0">的<span class="_ _6"> </span><span class="ff3">11*21<span class="_ _6"> </span></span>组<span class="_ _7"> </span>分别求出与之对应的<span class="_ _8"> </span>。</div><div class="t m0 x2 h5 y8 ff1 fs2 fc0 sc1 ls0 ws0">（<span class="ff3">2<span class="_ _1"></span></span>）<span class="_ _1"></span><span class="ff3">.</span>分<span class="_ _1"></span>片<span class="_ _1"></span>二<span class="_ _1"></span>次代<span class="_ _1"></span>数插<span class="_ _1"></span>值<span class="_ _1"></span>：对<span class="_ _1"></span>于<span class="_ _1"></span>已求<span class="_ _1"></span>出<span class="_ _1"></span>的<span class="_ _9"> </span>，<span class="_ _1"></span>使用<span class="_ _1"></span>分片<span class="_ _1"></span>二<span class="_ _1"></span>次<span class="_ _1"></span>代数<span class="_ _1"></span>插值<span class="_ _1"></span>法<span class="_ _1"></span>对原</div><div class="t m0 x2 h5 y9 ff1 fs2 fc0 sc1 ls0 ws0">题<span class="_ _a"></span>中<span class="_ _a"></span>的<span class="_ _b"> </span>的<span class="_ _a"></span>数<span class="_ _a"></span>表<span class="_ _a"></span>进<span class="_ _a"></span>行<span class="_ _a"></span>插<span class="_ _a"></span>值<span class="_ _a"></span>得<span class="_ _a"></span>到<span class="_ _a"></span>与<span class="_ _c"> </span>对<span class="_ _a"></span>应<span class="_ _a"></span>的<span class="_ _d"> </span>。<span class="_ _a"></span>于<span class="_ _a"></span>是<span class="_ _a"></span>产<span class="_ _a"></span>生<span class="_ _a"></span>了<span class="_ _e"> </span><span class="ff3">z=f(x,y)<span class="_ _a"></span></span>的</div><div class="t m0 x2 h5 ya ff3 fs2 fc0 sc1 ls0 ws0">11*21<span class="_ _6"> </span><span class="ff1">个数值解。</span></div><div class="t m0 x2 h5 yb ff1 fs2 fc0 sc1 ls0 ws0">（<span class="ff3">3</span>）<span class="ff3">.</span>二元曲面的最小二乘拟合：从<span class="_ _f"> </span><span class="ff3">k=1<span class="_ _6"> </span></span>开始逐渐增大<span class="_ _6"> </span><span class="ff3">k<span class="_ _6"> </span></span>的值，对<span class="_ _6"> </span><span class="ff3">z=f(x,y)</span>进</div><div class="t m0 x2 h5 yc ff1 fs2 fc0 sc1 ls0 ws0">行拟合，得到<span class="_ _1"></span>每次的<span class="_ _10"> </span>。当<span class="_ _11"> </span>时结束计算，<span class="_ _1"></span>输出拟合结果。（注：<span class="_ _1"></span>在计</div><div class="t m0 x2 h5 yd ff1 fs2 fc0 sc1 ls0 ws0">算中<span class="_ _1"></span>为了<span class="_ _1"></span>避免<span class="_ _1"></span>矩阵<span class="_ _1"></span>求<span class="_ _1"></span>逆，<span class="_ _1"></span>把在<span class="_ _1"></span>矩阵<span class="_ _1"></span>中需<span class="_ _1"></span>要求<span class="_ _1"></span>逆的<span class="_ _1"></span>部分<span class="_ _1"></span>转换<span class="_ _1"></span>成解<span class="_ _1"></span>线性<span class="_ _1"></span>方程<span class="_ _1"></span>组<span class="_ _1"></span>如：</div><div class="t m0 x4 h5 ye ff1 fs2 fc0 sc1 ls0 ws0">可<span class="_ _a"></span>以<span class="_ _a"></span>令<span class="_ _12"> </span>推<span class="_ _a"></span>得<span class="_ _13"> </span>再<span class="_ _a"></span>用<span class="_ _a"></span>同<span class="_ _a"></span>样<span class="_ _a"></span>方<span class="_ _a"></span>法</div><div class="t m0 x2 h5 yf ff1 fs2 fc0 sc1 ls0 ws0">计算<span class="_ _14"> </span>中<span class="_ _1"></span>的<span class="_ _6"> </span><span class="ff3">C<span class="_ _15"> </span>…<span class="_ _1"></span>…<span class="_ _16"></span><span class="ff1">（<span class="_ _1"></span>可以<span class="_ _1"></span>先转<span class="_ _1"></span>置<span class="_ _17"> </span>）<span class="_ _1"></span>，本<span class="_ _1"></span>程序<span class="_ _1"></span>中<span class="_ _1"></span>解<span class="_ _1"></span>线性<span class="_ _1"></span>方程<span class="_ _1"></span>组<span class="_ _1"></span>所采<span class="_ _1"></span>用</span></span></div><div class="t m0 x2 h5 y10 ff1 fs2 fc0 sc1 ls0 ws0">的方法为列主元高斯消元法。）</div><div class="t m0 x2 h5 y11 ff3 fs2 fc0 sc1 ls0 ws0">4.<span class="_ _1"></span><span class="ff1">观<span class="_ _1"></span>察<span class="_ _a"></span>逼<span class="_ _1"></span>近<span class="_ _1"></span>效<span class="_ _1"></span>果<span class="_ _a"></span>：<span class="_ _1"></span>计<span class="_ _1"></span>算<span class="_ _18"> </span>的<span class="_ _1"></span>值<span class="_ _1"></span>并<span class="_ _a"></span>输<span class="_ _1"></span>出<span class="_ _1"></span>结</span></div><div class="t m0 x2 h5 y12 ff1 fs2 fc0 sc1 ls0 ws0">果，以观察逼近的效果。其中<span class="_ _19"> </span>。</div><div class="t m0 x2 h4 y13 ff2 fs1 fc0 sc1 ls0 ws0">2.<span class="_ _0"> </span><span class="ff1">主要函数功能：</span></div><div class="t m0 x2 h6 y14 ff3 fs2 fc0 sc1 ls0 ws0">NewtonIter(double x, double y, vector< double ><span class="_ _1"></span> &x0)</div><div class="t m0 x2 h5 y15 ff1 fs2 fc0 sc1 ls0 ws0">：用来进行牛顿迭代法解线性方程；</div><div class="t m0 x2 h6 y16 ff3 fs2 fc0 sc1 ls0 ws0">SurfFit(vector< double > xi, vector< double > yj,</div><div class="t m0 x5 h6 y17 ff3 fs2 fc0 sc1 ls0 ws0"> vector< vector< double > > U, int k)</div><div class="t m0 x2 h5 y18 ff1 fs2 fc0 sc1 ls0 ws0">：曲面最小二乘拟合；</div><div class="t m0 x2 h6 y19 ff3 fs2 fc0 sc1 ls0 ws0">QuadratiInterp(double x, double y, vector<double> u,</div><div class="t m0 x6 h6 y1a ff3 fs2 fc0 sc1 ls0 ws0"> vector<double> t, vector<vector<double>> z)</div><div class="t m0 x2 h5 y1b ff1 fs2 fc0 sc1 ls0 ws0">：分片线性二次插值；</div><div class="t m0 x2 h5 y1c ff1 fs2 fc0 sc1 ls0 ws0">其他小函数的详细说明见代码注释。</div><div class="t m0 x2 h4 y1d ff3 fs1 fc0 sc1 ls0 ws0">3.<span class="_ _1a"> </span><span class="ff1">源代码</span></div><div class="t m0 x2 h5 y1e ff3 fs3 fc0 sc1 ls0 ws0"> <span class="fs2">// MifaTest.cpp : <span class="ff1">定义控制台应用程序的入口点。</span></span></div><div class="t m0 x2 h6 y1f ff3 fs2 fc0 sc1 ls0 ws0">//</div><div class="t m0 x2 h5 y20 ff3 fs2 fc0 sc1 ls0 ws0">// QuadraticFit.cpp : <span class="ff1">定义控制台应用程序的入口点。</span></div><div class="t m0 x2 h6 y21 ff3 fs2 fc0 sc1 ls0 ws0">//</div><div class="t m0 x2 h6 y22 ff3 fs2 fc0 sc1 ls0 ws0">/</div><div class="t m0 x2 h6 y23 ff3 fs2 fc0 sc1 ls0 ws0">*********************************************************************</div></div><div class="t m0 x7 h7 y24 ff2 fs4 fc0 sc1 ls0 ws0">1</div></div><div class="pi" data-data='{"ctm":[1.611850,0.000000,0.000000,1.611850,0.000000,0.000000]}'></div></div>
</body>
</html>
<div id="pf2" class="pf w0 h0" data-page-no="2"><div class="pc pc2 w0 h0"><img class="bi x0 y0 w1 h1" alt="" src="https://static.pudn.com/prod/directory_preview_static/6273a38a40256a40ce34026f/bg2.jpg"><div class="c x0 y1 w2 h2"><div class="t m0 x2 h6 y25 ff3 fs2 fc0 sc1 ls0 ws0">*</div><div class="t m0 x2 h6 y26 ff3 fs2 fc0 sc1 ls0 ws0">Author: </div><div class="t m0 x2 h6 y27 ff3 fs2 fc0 sc1 ls0 ws0">Function List: </div><div class="t m0 x2 h6 y5 ff3 fs2 fc0 sc1 ls0 ws0">double Finf(vector< double > vec)</div><div class="t m0 x2 h5 y28 ff3 fs2 fc0 sc1 ls0 ws0"> :<span class="ff1">求解无穷范数</span></div><div class="t m0 x2 h6 y29 ff3 fs2 fc0 sc1 ls0 ws0">vector< double > GausSolv(vector< vector< double > > a,</div><div class="t m0 x8 h6 y2a ff3 fs2 fc0 sc1 ls0 ws0"> vector< double > b)</div><div class="t m0 x2 h5 y2b ff3 fs2 fc0 sc1 ls0 ws0"> <span class="ff1">：在插值时用的</span>Gauss<span class="ff1">法求解线性方程组函数</span></div><div class="t m0 x2 h6 y2c ff3 fs2 fc0 sc1 ls0 ws0">vector< double> NewtonIter(double x, double y, ve<span class="_ _1"></span>ctor< double > </div><div class="t m0 x2 h6 y2d ff3 fs2 fc0 sc1 ls0 ws0">&x0)</div><div class="t m0 x2 h5 y2e ff3 fs2 fc0 sc1 ls0 ws0"> <span class="ff1">：连续牛顿迭代法求解非线性方程组</span></div><div class="t m0 x2 h6 y2f ff3 fs2 fc0 sc1 ls0 ws0">int PosPoint(double val, vector<double> p)</div><div class="t m0 x2 h5 ya ff3 fs2 fc0 sc1 ls0 ws0"> <span class="ff1">：得到在插值中</span>t<span class="ff1">，</span>u<span class="ff1">所对应的节点的位置</span></div><div class="t m0 x2 h6 yb ff3 fs2 fc0 sc1 ls0 ws0">double L(int s, double val, int point, <span class="_ _1"></span>vector<span class="_ _1b"></span><double> vec)</div><div class="t m0 x2 h5 y30 ff3 fs2 fc0 sc1 ls0 ws0"> <span class="ff1">：根据节点位置信息计算</span>Lagrange<span class="ff1">插值基函数</span></div><div class="t m0 x2 h6 y31 ff3 fs2 fc0 sc1 ls0 ws0">double QuadratiInterp(double x, double y, vector<double> u,</div><div class="t m0 x9 h6 yd ff3 fs2 fc0 sc1 ls0 ws0"> vector<double> t, vector<vector<double>> z)</div><div class="t m0 x2 h5 y32 ff3 fs2 fc0 sc1 ls0 ws0"> <span class="ff1">：分片二次插值</span></div><div class="t m0 x2 h6 y33 ff3 fs2 fc0 sc1 ls0 ws0">vector< vector< double > > MatrixT(vector< vector< double > ><span class="_ _1"></span> a)</div><div class="t m0 x2 h5 y34 ff3 fs2 fc0 sc1 ls0 ws0"> <span class="ff1">：矩阵的转置</span></div><div class="t m0 x2 h6 y35 ff3 fs2 fc0 sc1 ls0 ws0">vector< vector< double > > ConveArraMulitArra(vector< vector< </div><div class="t m0 x2 h6 y10 ff3 fs2 fc0 sc1 ls0 ws0">double > >a, </div><div class="t m0 x2 h6 y36 ff3 fs2 fc0 sc1 ls0 ws0"> <span class="_ _1"></span> vector<<span class="_ _1"></span> vector< double > > b)</div><div class="t m0 x2 h5 y37 ff3 fs2 fc0 sc1 ls0 ws0"> <span class="ff1">：矩阵转置乘另一个矩阵得到新的矩阵</span></div><div class="t m0 x2 h6 y38 ff3 fs2 fc0 sc1 ls0 ws0">vector< vector< double > > GausElim(vector< vector< double > > </div><div class="t m0 x2 h6 y39 ff3 fs2 fc0 sc1 ls0 ws0">a,</div><div class="t m0 xa h6 y3a ff3 fs2 fc0 sc1 ls0 ws0">vector< vector< double > > b)</div><div class="t m0 x2 h5 y3b ff3 fs2 fc0 sc1 ls0 ws0"> <span class="ff1">：列主元</span>Gauss<span class="ff1">消元法，专门针对矩阵乘未知矩阵得到一个矩阵（需要求解闻</span></div><div class="t m0 x2 h5 y14 ff1 fs2 fc0 sc1 ls0 ws0">之矩阵）</div><div class="t m0 x2 h5 y15 ff3 fs2 fc0 sc1 ls0 ws0"> <span class="ff1">这种方法比一次次调用普通的得到解向量再组合有更高的效率</span></div><div class="t m0 x2 h6 y16 ff3 fs2 fc0 sc1 ls0 ws0">vector< vector< double > > SurfFit(vector< double > xi, vector<<span class="_ _1"></span> </div><div class="t m0 x2 h6 y17 ff3 fs2 fc0 sc1 ls0 ws0">double > yj,</div><div class="t m0 xb h6 y18 ff3 fs2 fc0 sc1 ls0 ws0"> vector< vector< double > > U, int <span class="_ _1"></span>k)</div><div class="t m0 x2 h5 y19 ff3 fs2 fc0 sc1 ls0 ws0"> <span class="ff1">：曲面拟合</span></div><div class="t m0 x2 h6 y1a ff3 fs2 fc0 sc1 ls0 ws0">vector< vector< double > > GetPxyArra(vector< double > x, </div><div class="t m0 x2 h6 y1b ff3 fs2 fc0 sc1 ls0 ws0">vector< double > y, </div><div class="t m0 x2 h6 y1c ff3 fs2 fc0 sc1 ls0 ws0"> <span class="_ _1"></span> vector< vector< double ><span class="_ _1"></span> > C)</div><div class="t m0 x2 h5 y3c ff3 fs2 fc0 sc1 ls0 ws0"> <span class="ff1">：根据拟合出来的系数矩阵</span>C<span class="ff1">求解</span>p<span class="ff1">（</span>x<span class="ff1">，</span>y<span class="ff1">）</span></div><div class="t m0 x2 h6 y3d ff3 fs2 fc0 sc1 ls0 ws0">double CacuDeta(vector< vector< double > > f, vector< vector<<span class="_ _1"></span> </div><div class="t m0 x2 h6 y1e ff3 fs2 fc0 sc1 ls0 ws0">double > > p)</div><div class="t m0 x2 h5 y1f ff3 fs2 fc0 sc1 ls0 ws0"> <span class="ff1">：计算误差条件</span>σ<span class="ff1">的函数</span></div><div class="t m0 x2 h6 y21 ff3 fs2 fc0 sc1 ls0 ws0">*********************************************************************</div><div class="t m0 x2 h6 y22 ff3 fs2 fc0 sc1 ls0 ws0">*******/</div></div><div class="t m0 x7 h7 y24 ff2 fs4 fc0 sc1 ls0 ws0">2</div></div><div class="pi" data-data='{"ctm":[1.611850,0.000000,0.000000,1.611850,0.000000,0.000000]}'></div></div>
<div id="pf3" class="pf w0 h0" data-page-no="3"><div class="pc pc3 w0 h0"><img class="bi x0 y0 w1 h1" alt="" src="https://static.pudn.com/prod/directory_preview_static/6273a38a40256a40ce34026f/bg3.jpg"><div class="c x0 y1 w2 h2"><div class="t m0 x2 h6 y3e ff3 fs2 fc0 sc1 ls0 ws0">#include "stdafx.h"</div><div class="t m0 x2 h6 y26 ff3 fs2 fc0 sc1 ls0 ws0">#include <vector></div><div class="t m0 x2 h6 y27 ff3 fs2 fc0 sc1 ls0 ws0">#include <cmath></div><div class="t m0 x2 h6 y4 ff3 fs2 fc0 sc1 ls0 ws0">#include <iostream></div><div class="t m0 x2 h6 y5 ff3 fs2 fc0 sc1 ls0 ws0">#include <fstream></div><div class="t m0 x2 h6 y28 ff3 fs2 fc0 sc1 ls0 ws0">using namespace std;</div><div class="t m0 x2 h5 y2a ff3 fs2 fc0 sc1 ls0 ws0">/*****<span class="ff1">求无穷范数</span>********/</div><div class="t m0 x2 h6 y2b ff3 fs2 fc0 sc1 ls0 ws0">double Finf(vector< double > vec)</div><div class="t m0 x2 h6 y2c ff3 fs2 fc0 sc1 ls0 ws0">{</div><div class="t m0 xc h6 y2d ff3 fs2 fc0 sc1 ls0 ws0">int n = vec.size();</div><div class="t m0 xc h6 y2e ff3 fs2 fc0 sc1 ls0 ws0">double temp = fabs(vec.at(0));</div><div class="t m0 xc h6 y2f ff3 fs2 fc0 sc1 ls0 ws0">for (int i=1; i<n; i++)</div><div class="t m0 xc h6 ya ff3 fs2 fc0 sc1 ls0 ws0">{</div><div class="t m0 x5 h6 yb ff3 fs2 fc0 sc1 ls0 ws0">if (fabs(vec.at(i)) > temp)</div><div class="t m0 x5 h6 y30 ff3 fs2 fc0 sc1 ls0 ws0">{</div><div class="t m0 xd h6 y31 ff3 fs2 fc0 sc1 ls0 ws0">temp = fabs(vec.at(i));</div><div class="t m0 x5 h6 yd ff3 fs2 fc0 sc1 ls0 ws0">}</div><div class="t m0 xc h6 y32 ff3 fs2 fc0 sc1 ls0 ws0">}</div><div class="t m0 xc h6 y33 ff3 fs2 fc0 sc1 ls0 ws0">return temp;</div><div class="t m0 x2 h6 y34 ff3 fs2 fc0 sc1 ls0 ws0">}</div><div class="t m0 x2 h6 y35 ff3 fs2 fc0 sc1 ls0 ws0">vector< double > VectSubs(vector< double > &a, vector< double ></div><div class="t m0 x2 h6 y10 ff3 fs2 fc0 sc1 ls0 ws0">&b)</div><div class="t m0 x2 h6 y36 ff3 fs2 fc0 sc1 ls0 ws0">{</div><div class="t m0 xc h6 y37 ff3 fs2 fc0 sc1 ls0 ws0">int s1 = a.size();</div><div class="t m0 xc h6 y38 ff3 fs2 fc0 sc1 ls0 ws0">int s2 = b.size();</div><div class="t m0 xc h6 y39 ff3 fs2 fc0 sc1 ls0 ws0">if(s1 != s2)</div><div class="t m0 xc h6 y3a ff3 fs2 fc0 sc1 ls0 ws0">{</div><div class="t m0 x5 h6 y3b ff3 fs2 fc0 sc1 ls0 ws0">cerr << "Vectors not match in size !";</div><div class="t m0 xc h6 y14 ff3 fs2 fc0 sc1 ls0 ws0">}</div><div class="t m0 xc h6 y15 ff3 fs2 fc0 sc1 ls0 ws0">else</div><div class="t m0 xc h6 y16 ff3 fs2 fc0 sc1 ls0 ws0">{</div><div class="t m0 x5 h6 y17 ff3 fs2 fc0 sc1 ls0 ws0">vector< double > value(s1,0);</div><div class="t m0 x5 h6 y18 ff3 fs2 fc0 sc1 ls0 ws0">for (int i=1;i<s1;i++)</div><div class="t m0 x5 h6 y19 ff3 fs2 fc0 sc1 ls0 ws0">{</div><div class="t m0 xd h6 y1a ff3 fs2 fc0 sc1 ls0 ws0">value[i] = a[i] - b[i];</div><div class="t m0 x5 h6 y1b ff3 fs2 fc0 sc1 ls0 ws0">}</div><div class="t m0 x5 h6 y1c ff3 fs2 fc0 sc1 ls0 ws0">return value;</div><div class="t m0 xc h6 y3c ff3 fs2 fc0 sc1 ls0 ws0">}</div><div class="t m0 x2 h6 y3d ff3 fs2 fc0 sc1 ls0 ws0">}</div><div class="t m0 x2 h6 y20 ff3 fs2 fc0 sc1 ls0 ws0">vector< double > GausSolv(vector< vector< double > > a,</div><div class="t m0 x8 h6 y21 ff3 fs2 fc0 sc1 ls0 ws0"> vector< double > b)</div><div class="t m0 x2 h6 y22 ff3 fs2 fc0 sc1 ls0 ws0">{</div><div class="t m0 xc h6 y23 ff3 fs2 fc0 sc1 ls0 ws0">int am = a.size();</div></div><div class="t m0 x7 h7 y24 ff2 fs4 fc0 sc1 ls0 ws0">3</div></div><div class="pi" data-data='{"ctm":[1.611850,0.000000,0.000000,1.611850,0.000000,0.000000]}'></div></div>
<div id="pf4" class="pf w0 h0" data-page-no="4"><div class="pc pc4 w0 h0"><img class="bi x0 y0 w1 h1" alt="" src="https://static.pudn.com/prod/directory_preview_static/6273a38a40256a40ce34026f/bg4.jpg"><div class="c x0 y1 w2 h2"><div class="t m0 xc h6 y25 ff3 fs2 fc0 sc1 ls0 ws0">int an = a.at(0).size();</div><div class="t m0 xc h6 y3e ff3 fs2 fc0 sc1 ls0 ws0">int bm = b.size();</div><div class="t m0 xc h6 y26 ff3 fs2 fc0 sc1 ls0 ws0">if (am != bm)</div><div class="t m0 xc h6 y27 ff3 fs2 fc0 sc1 ls0 ws0">{</div><div class="t m0 x5 h6 y4 ff3 fs2 fc0 sc1 ls0 ws0">cerr<<"Array not match in size!!";</div><div class="t m0 xc h6 y5 ff3 fs2 fc0 sc1 ls0 ws0">}</div><div class="t m0 xc h6 y28 ff3 fs2 fc0 sc1 ls0 ws0">else</div><div class="t m0 xc h6 y29 ff3 fs2 fc0 sc1 ls0 ws0">{</div><div class="t m0 x5 h6 y2a ff3 fs2 fc0 sc1 ls0 ws0">vector< double > X(an,0);</div><div class="t m0 x5 h6 y2b ff3 fs2 fc0 sc1 ls0 ws0">for (int k = 0; k < am; k++)</div><div class="t m0 x5 h6 y2c ff3 fs2 fc0 sc1 ls0 ws0">{</div><div class="t m0 xd h6 y2d ff3 fs2 fc0 sc1 ls0 ws0">if (a[k][k]>=-1e-12&& a[k][k]<=1e-12)</div><div class="t m0 xd h6 y2e ff3 fs2 fc0 sc1 ls0 ws0">{</div><div class="t m0 x6 h5 y2f ff3 fs2 fc0 sc1 ls0 ws0">cerr << "<span class="ff1">矩阵无法求解！！</span>";</div><div class="t m0 x6 h6 ya ff3 fs2 fc0 sc1 ls0 ws0">exit(-1);</div><div class="t m0 xd h6 yb ff3 fs2 fc0 sc1 ls0 ws0">}</div><div class="t m0 xd h6 y30 ff3 fs2 fc0 sc1 ls0 ws0">for (int i=k+1; i<am; i++)</div><div class="t m0 xd h6 y31 ff3 fs2 fc0 sc1 ls0 ws0">{</div><div class="t m0 x6 h6 yd ff3 fs2 fc0 sc1 ls0 ws0">double m = a[i][k] / a[k][k];</div><div class="t m0 x6 h6 y32 ff3 fs2 fc0 sc1 ls0 ws0">for (int j=k+1; j<am; j++)</div><div class="t m0 x6 h6 y33 ff3 fs2 fc0 sc1 ls0 ws0">{</div><div class="t m0 x9 h6 y34 ff3 fs2 fc0 sc1 ls0 ws0">a[i][j] = a[i][j] - m * a[k][j];</div><div class="t m0 x6 h6 y35 ff3 fs2 fc0 sc1 ls0 ws0">}</div><div class="t m0 x6 h6 y36 ff3 fs2 fc0 sc1 ls0 ws0">b[i] = b[i] - m * b[k];</div><div class="t m0 xd h6 y38 ff3 fs2 fc0 sc1 ls0 ws0">}//end of for (int i=k+1; i<am...</div><div class="t m0 x5 h6 y39 ff3 fs2 fc0 sc1 ls0 ws0">}</div><div class="t m0 x5 h5 y3a ff3 fs2 fc0 sc1 ls0 ws0">//*********<span class="ff1">回代过程</span>*********</div><div class="t m0 x5 h6 y3b ff3 fs2 fc0 sc1 ls0 ws0">X[am-1] = b[am-1] / a[am-1][am-1];</div><div class="t m0 x5 h6 y14 ff3 fs2 fc0 sc1 ls0 ws0">for (int k=am-2; k>=0; k--)</div><div class="t m0 x5 h6 y15 ff3 fs2 fc0 sc1 ls0 ws0">{</div><div class="t m0 xd h6 y16 ff3 fs2 fc0 sc1 ls0 ws0">double temp = 0.0;</div><div class="t m0 xd h6 y17 ff3 fs2 fc0 sc1 ls0 ws0">for(int j=k+1; j<am; j++)</div><div class="t m0 xd h6 y18 ff3 fs2 fc0 sc1 ls0 ws0">{</div><div class="t m0 x6 h6 y19 ff3 fs2 fc0 sc1 ls0 ws0">temp += a[k][j] * X[j];</div><div class="t m0 xd h6 y1a ff3 fs2 fc0 sc1 ls0 ws0">}</div><div class="t m0 xd h6 y1b ff3 fs2 fc0 sc1 ls0 ws0">X[k] = (b[k] - temp) / a[k][k];</div><div class="t m0 x5 h6 y1c ff3 fs2 fc0 sc1 ls0 ws0">}</div><div class="t m0 x5 h6 y3d ff3 fs2 fc0 sc1 ls0 ws0">return X;</div><div class="t m0 xc h6 y1e ff3 fs2 fc0 sc1 ls0 ws0">}// end of else</div><div class="t m0 x2 h6 y20 ff3 fs2 fc0 sc1 ls0 ws0">}</div></div><div class="t m0 x7 h7 y24 ff2 fs4 fc0 sc1 ls0 ws0">4</div></div><div class="pi" data-data='{"ctm":[1.611850,0.000000,0.000000,1.611850,0.000000,0.000000]}'></div></div>
<div id="pf5" class="pf w0 h0" data-page-no="5"><div class="pc pc5 w0 h0"><img class="bi x0 y0 w1 h1" alt="" src="https://static.pudn.com/prod/directory_preview_static/6273a38a40256a40ce34026f/bg5.jpg"><div class="c x0 y1 w2 h2"><div class="t m0 x2 h6 y25 ff3 fs2 fc0 sc1 ls0 ws0">vector< double> NewtonIter(double x, double y, ve<span class="_ _1"></span>ctor< double > </div><div class="t m0 x2 h6 y3e ff3 fs2 fc0 sc1 ls0 ws0">&x0)</div><div class="t m0 x2 h6 y26 ff3 fs2 fc0 sc1 ls0 ws0">{</div><div class="t m0 xc h6 y27 ff3 fs2 fc0 sc1 ls0 ws0">vector< double > Fxk(4,0);</div><div class="t m0 xc h6 y4 ff3 fs2 fc0 sc1 ls0 ws0">vector< vector< double > > Dfxk(4,vector<double>(4,0));</div><div class="t m0 xc h6 y5 ff3 fs2 fc0 sc1 ls0 ws0">vector< double > deltax(4,0);</div><div class="t m0 xc h6 y28 ff3 fs2 fc0 sc1 ls0 ws0">vector< double > xk = x0;</div><div class="t m0 xc h6 y29 ff3 fs2 fc0 sc1 ls0 ws0">do</div><div class="t m0 xc h6 y2a ff3 fs2 fc0 sc1 ls0 ws0">{</div><div class="t m0 x5 h6 y2b ff3 fs2 fc0 sc1 ls0 ws0">for (int i=0; i<4; i++)</div><div class="t m0 x5 h6 y2c ff3 fs2 fc0 sc1 ls0 ws0">{</div><div class="t m0 xd h6 y2e ff3 fs2 fc0 sc1 ls0 ws0">xk.at(i) += deltax.at(i);</div><div class="t m0 x5 h6 ya ff3 fs2 fc0 sc1 ls0 ws0">}</div><div class="t m0 x5 h5 yb ff3 fs2 fc0 sc1 ls0 ws0">//<span class="ff1">计算</span>x k+1 = xk + delta x;</div><div class="t m0 x5 h5 y30 ff3 fs2 fc0 sc1 ls0 ws0">//<span class="ff1">计算</span>F<span class="ff1">（</span>xk<span class="ff1">）</span></div><div class="t m0 x5 h6 y31 ff3 fs2 fc0 sc1 ls0 ws0">Fxk.at(0) = -1.0*(0.5*cos(xk.at(1)) +xk.at(0) + xk.at(2) + </div><div class="t m0 x2 h6 yd ff3 fs2 fc0 sc1 ls0 ws0">xk.at(3) - x - 2.67);</div><div class="t m0 x5 h6 y32 ff3 fs2 fc0 sc1 ls0 ws0">Fxk.at(1) = -1.0*(xk.at(1) + 0.5*sin(xk.at(0)) + xk.at(2) + </div><div class="t m0 x2 h6 y33 ff3 fs2 fc0 sc1 ls0 ws0">xk.at(3) - y -1.07);</div><div class="t m0 x5 h6 y34 ff3 fs2 fc0 sc1 ls0 ws0">Fxk.at(2) = -1.0*(0.5*xk.at(1) + xk.at(0) + cos(xk.at(2)) + </div><div class="t m0 x2 h6 y35 ff3 fs2 fc0 sc1 ls0 ws0">xk.at(3) - x - 3.74);</div><div class="t m0 x5 h6 y10 ff3 fs2 fc0 sc1 ls0 ws0">Fxk.at(3) = -1.0*(xk.at(1) + 0.5*xk.at(0) + xk.at(2) + </div><div class="t m0 x2 h6 y36 ff3 fs2 fc0 sc1 ls0 ws0">sin(xk.at(3)) - y -0.79);</div><div class="t m0 x5 h5 y37 ff3 fs2 fc0 sc1 ls0 ws0">//<span class="ff1">给导数矩阵赋值</span></div><div class="t m0 x5 h6 y38 ff3 fs2 fc0 sc1 ls0 ws0">Dfxk[0][0] = 1;</div><div class="t m0 x5 h6 y39 ff3 fs2 fc0 sc1 ls0 ws0">Dfxk[0][1] = -0.5 * sin(xk.at(1));</div><div class="t m0 x5 h6 y3a ff3 fs2 fc0 sc1 ls0 ws0">Dfxk[0][2] = 1;</div><div class="t m0 x5 h6 y3b ff3 fs2 fc0 sc1 ls0 ws0">Dfxk[0][3] = 1;</div><div class="t m0 x5 h6 y15 ff3 fs2 fc0 sc1 ls0 ws0">Dfxk[1][0] = 0.5 * cos(xk.at(1));</div><div class="t m0 x5 h6 y16 ff3 fs2 fc0 sc1 ls0 ws0">Dfxk[1][1] = 1;</div><div class="t m0 x5 h6 y17 ff3 fs2 fc0 sc1 ls0 ws0">Dfxk[1][2] = 1;</div><div class="t m0 x5 h6 y18 ff3 fs2 fc0 sc1 ls0 ws0">Dfxk[1][3] = 1;</div><div class="t m0 x5 h6 y1a ff3 fs2 fc0 sc1 ls0 ws0">Dfxk[2][0] = 1;</div><div class="t m0 x5 h6 y1b ff3 fs2 fc0 sc1 ls0 ws0">Dfxk[2][1] = 0.5;</div><div class="t m0 x5 h6 y1c ff3 fs2 fc0 sc1 ls0 ws0">Dfxk[2][2] = -1.0*sin(xk.at(2));</div><div class="t m0 x5 h6 y3c ff3 fs2 fc0 sc1 ls0 ws0">Dfxk[2][3] = 1;</div><div class="t m0 x5 h6 y1e ff3 fs2 fc0 sc1 ls0 ws0">Dfxk[3][0] = 0.5;</div><div class="t m0 x5 h6 y1f ff3 fs2 fc0 sc1 ls0 ws0">Dfxk[3][1] = 1;</div><div class="t m0 x5 h6 y20 ff3 fs2 fc0 sc1 ls0 ws0">Dfxk[3][2] = 1;</div><div class="t m0 x5 h6 y21 ff3 fs2 fc0 sc1 ls0 ws0">Dfxk[3][3] = cos(xk.at(3));</div><div class="t m0 x5 h6 y23 ff3 fs2 fc0 sc1 ls0 ws0">deltax = GausSolv(Dfxk, Fxk);</div></div><div class="t m0 x7 h7 y24 ff2 fs4 fc0 sc1 ls0 ws0">5</div></div><div class="pi" data-data='{"ctm":[1.611850,0.000000,0.000000,1.611850,0.000000,0.000000]}'></div></div>
