Matlab一度被認為是最專業的數值計算工具之一,相信許多同學都或多或少用過這個工具。相比而言,Python作為一種膠水式的語言,其設計之初就不是為科學計算服務的。之前也看到許多人在吐槽說用Python去復現一些計算過程時經常失敗,因此(包括本人)也懷疑過是Python本身數值精度不夠導致的。那么Python的精度究竟如何,本文就來一探究竟。為了方便,我們就用線性方程組的求解來對比這一事實。
1、實驗設計
本文考慮的線性方程組: Ax=bAx=b A x = b
#mermaid-svg-KlQCVM1aHWdwKg9a .label {font-family: 'trebuchet ms', verdana, arial;font-family: var(--mermaid-font-family);fill: #333;color: #333; }#mermaid-svg-KlQCVM1aHWdwKg9a .label text {fill: #333; }#mermaid-svg-KlQCVM1aHWdwKg9a .node rect,
#mermaid-svg-KlQCVM1aHWdwKg9a .node circle,
#mermaid-svg-KlQCVM1aHWdwKg9a .node ellipse,
#mermaid-svg-KlQCVM1aHWdwKg9a .node polygon,
#mermaid-svg-KlQCVM1aHWdwKg9a .node path {fill: #ECECFF;stroke: #9370DB;stroke-width: 1px; }#mermaid-svg-KlQCVM1aHWdwKg9a .node .label {text-align: center;fill: #333; }#mermaid-svg-KlQCVM1aHWdwKg9a .node.clickable {cursor: pointer; }#mermaid-svg-KlQCVM1aHWdwKg9a .arrowheadPath {fill: #333333; }#mermaid-svg-KlQCVM1aHWdwKg9a .edgePath .path {stroke: #333333;stroke-width: 1.5px; }#mermaid-svg-KlQCVM1aHWdwKg9a .flowchart-link {stroke: #333333;fill: none; }#mermaid-svg-KlQCVM1aHWdwKg9a .edgeLabel {background-color: #e8e8e8;text-align: center; }
#mermaid-svg-KlQCVM1aHWdwKg9a .edgeLabel rect {opacity: 0.9; }
#mermaid-svg-KlQCVM1aHWdwKg9a .edgeLabel span {color: #333; }#mermaid-svg-KlQCVM1aHWdwKg9a .cluster rect {fill: #ffffde;stroke: #aaaa33;stroke-width: 1px; }#mermaid-svg-KlQCVM1aHWdwKg9a .cluster text {fill: #333; }#mermaid-svg-KlQCVM1aHWdwKg9a div.mermaidTooltip {position: absolute;text-align: center;max-width: 200px;padding: 2px;font-family: 'trebuchet ms', verdana, arial;font-family: var(--mermaid-font-family);font-size: 12px;background: #ffffde;border: 1px solid #aaaa33;border-radius: 2px;pointer-events: none;z-index: 100; }#mermaid-svg-KlQCVM1aHWdwKg9a .actor {stroke: #CCCCFF;fill: #ECECFF; }#mermaid-svg-KlQCVM1aHWdwKg9a text.actor > tspan {fill: black;stroke: none; }#mermaid-svg-KlQCVM1aHWdwKg9a .actor-line {stroke: grey; }#mermaid-svg-KlQCVM1aHWdwKg9a .messageLine0 {stroke-width: 1.5;stroke-dasharray: none;stroke: #333; }#mermaid-svg-KlQCVM1aHWdwKg9a .messageLine1 {stroke-width: 1.5;stroke-dasharray: 2, 2;stroke: #333; }#mermaid-svg-KlQCVM1aHWdwKg9a #arrowhead path {fill: #333;stroke: #333; }#mermaid-svg-KlQCVM1aHWdwKg9a .sequenceNumber {fill: white; }#mermaid-svg-KlQCVM1aHWdwKg9a #sequencenumber {fill: #333; }#mermaid-svg-KlQCVM1aHWdwKg9a #crosshead path {fill: #333;stroke: #333; }#mermaid-svg-KlQCVM1aHWdwKg9a .messageText {fill: #333;stroke: #333; }#mermaid-svg-KlQCVM1aHWdwKg9a .labelBox {stroke: #CCCCFF;fill: #ECECFF; }#mermaid-svg-KlQCVM1aHWdwKg9a .labelText,#mermaid-svg-KlQCVM1aHWdwKg9a .labelText > tspan {fill: black;stroke: none; }#mermaid-svg-KlQCVM1aHWdwKg9a .loopText,#mermaid-svg-KlQCVM1aHWdwKg9a .loopText > tspan {fill: black;stroke: none; }#mermaid-svg-KlQCVM1aHWdwKg9a .loopLine {stroke-width: 2px;stroke-dasharray: 2, 2;stroke: #CCCCFF;fill: #CCCCFF; }#mermaid-svg-KlQCVM1aHWdwKg9a .note {stroke: #aaaa33;fill: #fff5ad; }#mermaid-svg-KlQCVM1aHWdwKg9a .noteText,#mermaid-svg-KlQCVM1aHWdwKg9a .noteText > tspan {fill: black;stroke: none; }#mermaid-svg-KlQCVM1aHWdwKg9a .activation0 {fill: #f4f4f4;stroke: #666; }#mermaid-svg-KlQCVM1aHWdwKg9a .activation1 {fill: #f4f4f4;stroke: #666; }#mermaid-svg-KlQCVM1aHWdwKg9a .activation2 {fill: #f4f4f4;stroke: #666; }#mermaid-svg-KlQCVM1aHWdwKg9a .mermaid-main-font {font-family: "trebuchet ms", verdana, arial;font-family: var(--mermaid-font-family); }#mermaid-svg-KlQCVM1aHWdwKg9a .section {stroke: none;opacity: 0.2; }#mermaid-svg-KlQCVM1aHWdwKg9a .section0 {fill: rgba(102, 102, 255, 0.49); }#mermaid-svg-KlQCVM1aHWdwKg9a .section2 {fill: #fff400; }#mermaid-svg-KlQCVM1aHWdwKg9a .section1,
#mermaid-svg-KlQCVM1aHWdwKg9a .section3 {fill: white;opacity: 0.2; }#mermaid-svg-KlQCVM1aHWdwKg9a .sectionTitle0 {fill: #333; }#mermaid-svg-KlQCVM1aHWdwKg9a .sectionTitle1 {fill: #333; }#mermaid-svg-KlQCVM1aHWdwKg9a .sectionTitle2 {fill: #333; }#mermaid-svg-KlQCVM1aHWdwKg9a .sectionTitle3 {fill: #333; }#mermaid-svg-KlQCVM1aHWdwKg9a .sectionTitle {text-anchor: start;font-size: 11px;text-height: 14px;font-family: 'trebuchet ms', verdana, arial;font-family: var(--mermaid-font-family); }#mermaid-svg-KlQCVM1aHWdwKg9a .grid .tick {stroke: lightgrey;opacity: 0.8;shape-rendering: crispEdges; }
#mermaid-svg-KlQCVM1aHWdwKg9a .grid .tick text {font-family: 'trebuchet ms', verdana, arial;font-family: var(--mermaid-font-family); }#mermaid-svg-KlQCVM1aHWdwKg9a .grid path {stroke-width: 0; }#mermaid-svg-KlQCVM1aHWdwKg9a .today {fill: none;stroke: red;stroke-width: 2px; }#mermaid-svg-KlQCVM1aHWdwKg9a .task {stroke-width: 2; }#mermaid-svg-KlQCVM1aHWdwKg9a .taskText {text-anchor: middle;font-family: 'trebuchet ms', verdana, arial;font-family: var(--mermaid-font-family); }#mermaid-svg-KlQCVM1aHWdwKg9a .taskText:not([font-size]) {font-size: 11px; }#mermaid-svg-KlQCVM1aHWdwKg9a .taskTextOutsideRight {fill: black;text-anchor: start;font-size: 11px;font-family: 'trebuchet ms', verdana, arial;font-family: var(--mermaid-font-family); }#mermaid-svg-KlQCVM1aHWdwKg9a .taskTextOutsideLeft {fill: black;text-anchor: end;font-size: 11px; }#mermaid-svg-KlQCVM1aHWdwKg9a .task.clickable {cursor: pointer; }#mermaid-svg-KlQCVM1aHWdwKg9a .taskText.clickable {cursor: pointer;fill: #003163 !important;font-weight: bold; }#mermaid-svg-KlQCVM1aHWdwKg9a .taskTextOutsideLeft.clickable {cursor: pointer;fill: #003163 !important;font-weight: bold; }#mermaid-svg-KlQCVM1aHWdwKg9a .taskTextOutsideRight.clickable {cursor: pointer;fill: #003163 !important;font-weight: bold; }#mermaid-svg-KlQCVM1aHWdwKg9a .taskText0,
#mermaid-svg-KlQCVM1aHWdwKg9a .taskText1,
#mermaid-svg-KlQCVM1aHWdwKg9a .taskText2,
#mermaid-svg-KlQCVM1aHWdwKg9a .taskText3 {fill: white; }#mermaid-svg-KlQCVM1aHWdwKg9a .task0,
#mermaid-svg-KlQCVM1aHWdwKg9a .task1,
#mermaid-svg-KlQCVM1aHWdwKg9a .task2,
#mermaid-svg-KlQCVM1aHWdwKg9a .task3 {fill: #8a90dd;stroke: #534fbc; }#mermaid-svg-KlQCVM1aHWdwKg9a .taskTextOutside0,
#mermaid-svg-KlQCVM1aHWdwKg9a .taskTextOutside2 {fill: black; }#mermaid-svg-KlQCVM1aHWdwKg9a .taskTextOutside1,
#mermaid-svg-KlQCVM1aHWdwKg9a .taskTextOutside3 {fill: black; }#mermaid-svg-KlQCVM1aHWdwKg9a .active0,
#mermaid-svg-KlQCVM1aHWdwKg9a .active1,
#mermaid-svg-KlQCVM1aHWdwKg9a .active2,
#mermaid-svg-KlQCVM1aHWdwKg9a .active3 {fill: #bfc7ff;stroke: #534fbc; }#mermaid-svg-KlQCVM1aHWdwKg9a .activeText0,
#mermaid-svg-KlQCVM1aHWdwKg9a .activeText1,
#mermaid-svg-KlQCVM1aHWdwKg9a .activeText2,
#mermaid-svg-KlQCVM1aHWdwKg9a .activeText3 {fill: black !important; }#mermaid-svg-KlQCVM1aHWdwKg9a .done0,
#mermaid-svg-KlQCVM1aHWdwKg9a .done1,
#mermaid-svg-KlQCVM1aHWdwKg9a .done2,
#mermaid-svg-KlQCVM1aHWdwKg9a .done3 {stroke: grey;fill: lightgrey;stroke-width: 2; }#mermaid-svg-KlQCVM1aHWdwKg9a .doneText0,
#mermaid-svg-KlQCVM1aHWdwKg9a .doneText1,
#mermaid-svg-KlQCVM1aHWdwKg9a .doneText2,
#mermaid-svg-KlQCVM1aHWdwKg9a .doneText3 {fill: black !important; }#mermaid-svg-KlQCVM1aHWdwKg9a .crit0,
#mermaid-svg-KlQCVM1aHWdwKg9a .crit1,
#mermaid-svg-KlQCVM1aHWdwKg9a .crit2,
#mermaid-svg-KlQCVM1aHWdwKg9a .crit3 {stroke: #ff8888;fill: red;stroke-width: 2; }#mermaid-svg-KlQCVM1aHWdwKg9a .activeCrit0,
#mermaid-svg-KlQCVM1aHWdwKg9a .activeCrit1,
#mermaid-svg-KlQCVM1aHWdwKg9a .activeCrit2,
#mermaid-svg-KlQCVM1aHWdwKg9a .activeCrit3 {stroke: #ff8888;fill: #bfc7ff;stroke-width: 2; }#mermaid-svg-KlQCVM1aHWdwKg9a .doneCrit0,
#mermaid-svg-KlQCVM1aHWdwKg9a .doneCrit1,
#mermaid-svg-KlQCVM1aHWdwKg9a .doneCrit2,
#mermaid-svg-KlQCVM1aHWdwKg9a .doneCrit3 {stroke: #ff8888;fill: lightgrey;stroke-width: 2;cursor: pointer;shape-rendering: crispEdges; }#mermaid-svg-KlQCVM1aHWdwKg9a .milestone {transform: rotate(45deg) scale(0.8, 0.8); }#mermaid-svg-KlQCVM1aHWdwKg9a .milestoneText {font-style: italic; }#mermaid-svg-KlQCVM1aHWdwKg9a .doneCritText0,
#mermaid-svg-KlQCVM1aHWdwKg9a .doneCritText1,
#mermaid-svg-KlQCVM1aHWdwKg9a .doneCritText2,
#mermaid-svg-KlQCVM1aHWdwKg9a .doneCritText3 {fill: black !important; }#mermaid-svg-KlQCVM1aHWdwKg9a .activeCritText0,
#mermaid-svg-KlQCVM1aHWdwKg9a .activeCritText1,
#mermaid-svg-KlQCVM1aHWdwKg9a .activeCritText2,
#mermaid-svg-KlQCVM1aHWdwKg9a .activeCritText3 {fill: black !important; }#mermaid-svg-KlQCVM1aHWdwKg9a .titleText {text-anchor: middle;font-size: 18px;fill: black;font-family: 'trebuchet ms', verdana, arial;font-family: var(--mermaid-font-family); }#mermaid-svg-KlQCVM1aHWdwKg9a g.classGroup text {fill: #9370DB;stroke: none;font-family: 'trebuchet ms', verdana, arial;font-family: var(--mermaid-font-family);font-size: 10px; }
#mermaid-svg-KlQCVM1aHWdwKg9a g.classGroup text .title {font-weight: bolder; }#mermaid-svg-KlQCVM1aHWdwKg9a g.clickable {cursor: pointer; }#mermaid-svg-KlQCVM1aHWdwKg9a g.classGroup rect {fill: #ECECFF;stroke: #9370DB; }#mermaid-svg-KlQCVM1aHWdwKg9a g.classGroup line {stroke: #9370DB;stroke-width: 1; }#mermaid-svg-KlQCVM1aHWdwKg9a .classLabel .box {stroke: none;stroke-width: 0;fill: #ECECFF;opacity: 0.5; }#mermaid-svg-KlQCVM1aHWdwKg9a .classLabel .label {fill: #9370DB;font-size: 10px; }#mermaid-svg-KlQCVM1aHWdwKg9a .relation {stroke: #9370DB;stroke-width: 1;fill: none; }#mermaid-svg-KlQCVM1aHWdwKg9a .dashed-line {stroke-dasharray: 3; }#mermaid-svg-KlQCVM1aHWdwKg9a #compositionStart {fill: #9370DB;stroke: #9370DB;stroke-width: 1; }#mermaid-svg-KlQCVM1aHWdwKg9a #compositionEnd {fill: #9370DB;stroke: #9370DB;stroke-width: 1; }#mermaid-svg-KlQCVM1aHWdwKg9a #aggregationStart {fill: #ECECFF;stroke: #9370DB;stroke-width: 1; }#mermaid-svg-KlQCVM1aHWdwKg9a #aggregationEnd {fill: #ECECFF;stroke: #9370DB;stroke-width: 1; }#mermaid-svg-KlQCVM1aHWdwKg9a #dependencyStart {fill: #9370DB;stroke: #9370DB;stroke-width: 1; }#mermaid-svg-KlQCVM1aHWdwKg9a #dependencyEnd {fill: #9370DB;stroke: #9370DB;stroke-width: 1; }#mermaid-svg-KlQCVM1aHWdwKg9a #extensionStart {fill: #9370DB;stroke: #9370DB;stroke-width: 1; }#mermaid-svg-KlQCVM1aHWdwKg9a #extensionEnd {fill: #9370DB;stroke: #9370DB;stroke-width: 1; }#mermaid-svg-KlQCVM1aHWdwKg9a .commit-id,
#mermaid-svg-KlQCVM1aHWdwKg9a .commit-msg,
#mermaid-svg-KlQCVM1aHWdwKg9a .branch-label {fill: lightgrey;color: lightgrey;font-family: 'trebuchet ms', verdana, arial;font-family: var(--mermaid-font-family); }#mermaid-svg-KlQCVM1aHWdwKg9a .pieTitleText {text-anchor: middle;font-size: 25px;fill: black;font-family: 'trebuchet ms', verdana, arial;font-family: var(--mermaid-font-family); }#mermaid-svg-KlQCVM1aHWdwKg9a .slice {font-family: 'trebuchet ms', verdana, arial;font-family: var(--mermaid-font-family); }#mermaid-svg-KlQCVM1aHWdwKg9a g.stateGroup text {fill: #9370DB;stroke: none;font-size: 10px;font-family: 'trebuchet ms', verdana, arial;font-family: var(--mermaid-font-family); }#mermaid-svg-KlQCVM1aHWdwKg9a g.stateGroup text {fill: #9370DB;fill: #333;stroke: none;font-size: 10px; }#mermaid-svg-KlQCVM1aHWdwKg9a g.statediagram-cluster .cluster-label text {fill: #333; }#mermaid-svg-KlQCVM1aHWdwKg9a g.stateGroup .state-title {font-weight: bolder;fill: black; }#mermaid-svg-KlQCVM1aHWdwKg9a g.stateGroup rect {fill: #ECECFF;stroke: #9370DB; }#mermaid-svg-KlQCVM1aHWdwKg9a g.stateGroup line {stroke: #9370DB;stroke-width: 1; }#mermaid-svg-KlQCVM1aHWdwKg9a .transition {stroke: #9370DB;stroke-width: 1;fill: none; }#mermaid-svg-KlQCVM1aHWdwKg9a .stateGroup .composit {fill: white;border-bottom: 1px; }#mermaid-svg-KlQCVM1aHWdwKg9a .stateGroup .alt-composit {fill: #e0e0e0;border-bottom: 1px; }#mermaid-svg-KlQCVM1aHWdwKg9a .state-note {stroke: #aaaa33;fill: #fff5ad; }
#mermaid-svg-KlQCVM1aHWdwKg9a .state-note text {fill: black;stroke: none;font-size: 10px; }#mermaid-svg-KlQCVM1aHWdwKg9a .stateLabel .box {stroke: none;stroke-width: 0;fill: #ECECFF;opacity: 0.7; }#mermaid-svg-KlQCVM1aHWdwKg9a .edgeLabel text {fill: #333; }#mermaid-svg-KlQCVM1aHWdwKg9a .stateLabel text {fill: black;font-size: 10px;font-weight: bold;font-family: 'trebuchet ms', verdana, arial;font-family: var(--mermaid-font-family); }#mermaid-svg-KlQCVM1aHWdwKg9a .node circle.state-start {fill: black;stroke: black; }#mermaid-svg-KlQCVM1aHWdwKg9a .node circle.state-end {fill: black;stroke: white;stroke-width: 1.5; }#mermaid-svg-KlQCVM1aHWdwKg9a #statediagram-barbEnd {fill: #9370DB; }#mermaid-svg-KlQCVM1aHWdwKg9a .statediagram-cluster rect {fill: #ECECFF;stroke: #9370DB;stroke-width: 1px; }#mermaid-svg-KlQCVM1aHWdwKg9a .statediagram-cluster rect.outer {rx: 5px;ry: 5px; }#mermaid-svg-KlQCVM1aHWdwKg9a .statediagram-state .divider {stroke: #9370DB; }#mermaid-svg-KlQCVM1aHWdwKg9a .statediagram-state .title-state {rx: 5px;ry: 5px; }#mermaid-svg-KlQCVM1aHWdwKg9a .statediagram-cluster.statediagram-cluster .inner {fill: white; }#mermaid-svg-KlQCVM1aHWdwKg9a .statediagram-cluster.statediagram-cluster-alt .inner {fill: #e0e0e0; }#mermaid-svg-KlQCVM1aHWdwKg9a .statediagram-cluster .inner {rx: 0;ry: 0; }#mermaid-svg-KlQCVM1aHWdwKg9a .statediagram-state rect.basic {rx: 5px;ry: 5px; }#mermaid-svg-KlQCVM1aHWdwKg9a .statediagram-state rect.divider {stroke-dasharray: 10,10;fill: #efefef; }#mermaid-svg-KlQCVM1aHWdwKg9a .note-edge {stroke-dasharray: 5; }#mermaid-svg-KlQCVM1aHWdwKg9a .statediagram-note rect {fill: #fff5ad;stroke: #aaaa33;stroke-width: 1px;rx: 0;ry: 0; }:root {--mermaid-font-family: '"trebuchet ms", verdana, arial';--mermaid-font-family: "Comic Sans MS", "Comic Sans", cursive; }#mermaid-svg-KlQCVM1aHWdwKg9a .error-icon {fill: #552222; }#mermaid-svg-KlQCVM1aHWdwKg9a .error-text {fill: #552222;stroke: #552222; }#mermaid-svg-KlQCVM1aHWdwKg9a .edge-thickness-normal {stroke-width: 2px; }#mermaid-svg-KlQCVM1aHWdwKg9a .edge-thickness-thick {stroke-width: 3.5px; }#mermaid-svg-KlQCVM1aHWdwKg9a .edge-pattern-solid {stroke-dasharray: 0; }#mermaid-svg-KlQCVM1aHWdwKg9a .edge-pattern-dashed {stroke-dasharray: 3; }#mermaid-svg-KlQCVM1aHWdwKg9a .edge-pattern-dotted {stroke-dasharray: 2; }#mermaid-svg-KlQCVM1aHWdwKg9a .marker {fill: #333333; }#mermaid-svg-KlQCVM1aHWdwKg9a .marker.cross {stroke: #333333; }:root { --mermaid-font-family: "trebuchet ms", verdana, arial;}#mermaid-svg-KlQCVM1aHWdwKg9a {color: rgba(0, 0, 0, 0.75);font: ;}
生成隨機系數矩陣A和常數項b對數據進行截斷使得兩種程序中數據完全相等分別用Matlab和Python求解計算各自誤差值
這里第二步其實很重要,由于Matlab默認數據類型為double,而Python的默認類型為float,如果直接粘貼或者導入就有可能在Python中被截斷,從而導致在一開始使用的數據就存在誤差,使得實驗結果不可信。
第四步主要是計算兩種程序跑完后均方根誤差: ∥Ax?b∥\|Ax-b\| ∥ A x ? b ∥ 這種驗證方式比較簡單。之所以這么做是因為數值計算過程中不一定能夠得到100%的精度,這個問題的主要原因是在數值計算時數字不會一直保持為真實值。尤其在除法時由于數值長度的限制會被強行截斷一部分,因此用這個指標來驗證方法是否正確是很有必要的。當然,根據它的定義不難看出,這個值最小精度就越高。
2、數據生成和處理
首先由Matlab隨機生成一個矩陣和一個向量:
A = rand(5,5)
b = rand(5,1)
然后在Matlab中將數據粘至txt文件,再整理成Python的格式粘生成Numpy數組,同時也將數據重新粘回Matlab代碼。這里原始數據是這個樣子:
A的值
0.814723686393179 0.0975404049994095 0.157613081677548 0.141886338627215 0.655740699156587
0.905791937075619 0.278498218867048 0.970592781760616 0.421761282626275 0.0357116785741896
0.126986816293506 0.546881519204984 0.957166948242946 0.915735525189067 0.849129305868777
0.913375856139019 0.957506835434298 0.485375648722841 0.792207329559554 0.933993247757551
0.632359246225410 0.964888535199277 0.800280468888800 0.959492426392903 0.678735154857774
b的值
0.757700000000000
0.743100000000000
0.392200000000000
0.655500000000000
0.171200000000000
當然也可以采用round函數一類的方式進行處理。
注意 這種操作雖然看似簡單并且有點笨拙,但能很好地保證二者在初始化時得到的結果是一樣的。
Matlab的操作相對簡單,直接將txt中的內容首尾加上[],再粘回代碼即可。對于Python要稍加處理,最終的代碼是:
import numpy
as np
A
= np
. array
( [ 0.814723686393179 , 0.905791937075619 , 0.126986816293506 , 0.913375856139019 , 0.632359246225410 , 0.0975404049994095 , 0.278498218867048 , 0.546881519204984 , 0.957506835434298 , 0.964888535199277 , 0.157613081677548 , 0.970592781760616 , 0.957166948242946 , 0.485375648722841 , 0.800280468888800 , 0.141886338627215 , 0.421761282626275 , 0.915735525189067 , 0.792207329559554 , 0.959492426392903 , 0.655740699156587 , 0.0357116785741896 , 0.849129305868777 , 0.933993247757551 , 0.678735154857774 ] ) . reshape
( [ 5 , 5 ] ) . Tb
= np
. array
( [ 0.757700000000000 , 0.743100000000000 , 0.392200000000000 , 0.655500000000000 , 0.171200000000000 ] )
注意這里Python的代碼最后要進行一下轉置,否則與Matlab的數據就就不一致。
整理完成看看數據和txt中的形狀就完全一致了:
Aarray
( [ [ 0.81472369 , 0.0975404 , 0.15761308 , 0.14188634 , 0.6557407 ] , [ 0.90579194 , 0.27849822 , 0.97059278 , 0.42176128 , 0.03571168 ] , [ 0.12698682 , 0.54688152 , 0.95716695 , 0.91573553 , 0.84912931 ] , [ 0.91337586 , 0.95750684 , 0.48537565 , 0.79220733 , 0.93399325 ] , [ 0.63235925 , 0.96488854 , 0.80028047 , 0.95949243 , 0.67873515 ] ] ) barray
( [ 0.7577 , 0.7431 , 0.3922 , 0.6555 , 0.1712 ] )
這里b變成了一維數組但不影響計算。
3、求解方程組
用Matlab求解線性方程組十分簡單,直接一個反除號就可以搞定:
>> x = A \ bx =-0.84093.77013.8572-8.00502.4445
注意這里x跑出來的值只是顯示值,實際值的位數是遠大于這個值的:
-0.840928749143009 3.77006972003583 3.85716025378275 -8.00495386938441 2.44448015270832
用Python求解線性方程組的方式其實也很多,最簡單的辦法就是直接用scipy庫的linalg模塊中的solve方法:
from scipy
. linalg
import solvex
= solve
( A
, b
) xarray
( [ - 0.84092875 , 3.77006972 , 3.85716025 , - 8.00495387 , 2.44448015 ] )
為了更準確地對比數值的精度,我們將x的值用np.savetxt函數保存至文件再看具體數字:
- 8.409287491430091910e-01
3.770069720035832184e+00
3.857160253782746739e+00
- 8.004953869384410226e+00
2.444480152708315313e+00
仔細觀察發現,Python輸出的數值位數要稍多一些,這是因為我們用Matlab粘貼的時候它默認顯示的位數不太一樣。
注:這也是為什么我們一定要在最開始進行數據統一規范的原因。Matlab在存儲變量時數據是保存了完整的位數的,但即便從變量查看窗口粘出來的數值也不一定是完整的數字。
4、精度驗證
error = norm(A*x-b)error =1.5801e-15% 查看15位
% 1.580118849929844e-15
np
. sqrt
( np
. sum ( ( A
. dot
( x
) - b
) ** 2 ) ) 1.580118849929844e-15
注意,由于兩個數的量級都是1e-15,因此這里只需要看它們各自前16位即可。仔細觀察不難發現,兩個程序得到的結果完全一致。
這里我們把有效數字從第一位到最后一位排列直接放出來對比
1580118849929843973242195223453751011511480618378089335607228349545039236545562744140625
% matlab
1580118849929843973242195223453751011511480618378089335607228349545039236545562744140625
肉眼觀察信不過,全部放入Python字符串試試對不對:
s1
= '1580118849929843973242195223453751011511480618378089335607228349545039236545562744140625'
s2
= '1580118849929843973242195223453751011511480618378089335607228349545039236545562744140625'
判斷一下是否相等:
s1
== s2
True
搞定!
5、結論
Python(Numpy, Scipy)的數值精度與Matlab是完全一致的。
其本質的原因來自于二者在矩陣計算都是使用的LAPACK, 詳情可見: Matlab和Python(Numpy,Scipy)與Lapack的關系
總結
以上是生活随笔 為你收集整理的深度对比Python(Numpy,Scipy)与Matlab的数值精度 的全部內容,希望文章能夠幫你解決所遇到的問題。
如果覺得生活随笔 網站內容還不錯,歡迎將生活随笔 推薦給好友。