Add Xianjin Chen's computational methods code

syx777 · Jun 13, 2020 · 3add6e4 · 3add6e4
1 parent 4aef83c
commit 3add6e4
Show file tree

Hide file tree

Showing 35 changed files with 1,289 additions and 0 deletions.
diff --git a/计算方法/labs/1/1.cpp b/计算方法/labs/1/1.cpp
@@ -0,0 +1,53 @@
+#include <bits/stdc++.h>
+using namespace std;
+float f(float x)
+{
+	return sqrt(x * x + 4) - 2;
+}
+float g(float x)
+{
+	return x * x / (sqrt(x * x + 4) + 2);
+}
+
+int main(){
+	// Problem 1
+	float x;
+	x = 1.;
+	for(int i = 0; i < 10; i++) {
+		x = x / 8;
+		cout << setiosflags(ios::scientific) << setprecision(12) << x << "\t" << f(x) << "\t" << g(x) << endl;
+	}
+	// Problem 2
+	double arr[] = {4040.045551380452, -2759471.276702747, -31.64291531266504,  2755462.874010974,  0.0000557052996742893};
+	double res;
+	// seq
+	res = 0.;
+	for(int i = 0; i < 5; i++) {
+		res += arr[i];
+	}
+	cout << setprecision(7) << "Final " << res << endl;
+	// rev
+	res = 0.;
+	for(int i = 4; i >= 0; i--) {
+		res += arr[i];
+	}
+	cout << "Final " << res << endl;
+	// posi and nega
+	res = 0.;
+	res += arr[0];
+	res += arr[3];
+	res += arr[4];
+	res += arr[1];
+	res += arr[2];
+	cout << "Final " << res << endl;
+	// from small to big
+	res = 0.;
+	res += arr[4];
+	res += arr[2];
+	res += arr[0];
+	res += arr[3];
+	res += arr[1];
+	cout << "Final " << res << endl;
+	// Real result is 0.
+	return 0;
+}
diff --git a/计算方法/labs/1/1.docx b/计算方法/labs/1/1.docx
diff --git a/计算方法/labs/1/cm.nb b/计算方法/labs/1/cm.nb
@@ -0,0 +1,191 @@
+(* Content-type: application/vnd.wolfram.mathematica *)
+
+(*** Wolfram Notebook File ***)
+(* http://www.wolfram.com/nb *)
+
+(* CreatedBy='Mathematica 11.3' *)
+
+(*CacheID: 234*)
+(* Internal cache information:
+NotebookFileLineBreakTest
+NotebookFileLineBreakTest
+NotebookDataPosition[       158,          7]
+NotebookDataLength[      6052,        183]
+NotebookOptionsPosition[      5168,        158]
+NotebookOutlinePosition[      5533,        174]
+CellTagsIndexPosition[      5490,        171]
+WindowFrame->Normal*)
+
+(* Beginning of Notebook Content *)
+Notebook[{
+
+Cell[CellGroupData[{
+Cell[BoxData[
+ RowBox[{
+  RowBox[{"X", "=", 
+   RowBox[{"{", 
+    RowBox[{"4040.045551380452", ",", 
+     RowBox[{"-", "2759471.276702747"}], ",", 
+     RowBox[{"-", "31.64291531266504"}], ",", "2755462.874010974", ",", 
+     "0.0000557052996742893"}], "}"}]}], "\n"}]], "Input",
+ CellChangeTimes->{{3.7602515756829224`*^9, 3.760251609278879*^9}, {
+  3.7602516900422325`*^9, 3.7602516955699854`*^9}},
+ CellLabel->"In[8]:=",ExpressionUUID->"27aa795a-19b1-4daa-99c2-798100ecb41b"],
+
+Cell[BoxData[
+ RowBox[{"{", 
+  RowBox[{"4040.045551380452`", ",", 
+   RowBox[{"-", "2.759471276702747`*^6"}], ",", 
+   RowBox[{"-", "31.64291531266504`"}], ",", "2.755462874010974`*^6", ",", 
+   "0.0000557052996742893`"}], "}"}]], "Output",
+ CellChangeTimes->{
+  3.760251609675742*^9, {3.7602516904743037`*^9, 3.7602516960183907`*^9}},
+ CellLabel->"Out[8]=",ExpressionUUID->"3d644f9c-94d5-4c0c-8f97-2c028531f19d"]
+}, Open  ]],
+
+Cell[CellGroupData[{
+
+Cell[BoxData[
+ RowBox[{"Total", "[", "X", "]"}]], "Input",
+ CellChangeTimes->{{3.7602516142985*^9, 3.760251640514464*^9}},
+ CellLabel->"In[9]:=",ExpressionUUID->"abc47792-61a2-410c-851a-59cf59a6a8bd"],
+
+Cell[BoxData["0.`"], "Output",
+ CellChangeTimes->{{3.7602516171608987`*^9, 3.760251640820704*^9}, {
+  3.7602516932206087`*^9, 3.760251698447276*^9}},
+ CellLabel->"Out[9]=",ExpressionUUID->"6c2dff50-bcc1-49d8-a108-36311f38843c"]
+}, Open  ]],
+
+Cell[CellGroupData[{
+
+Cell[BoxData[
+ RowBox[{"Table", "[", 
+  RowBox[{
+   RowBox[{
+    RowBox[{"N", "[", 
+     RowBox[{
+      RowBox[{
+       SqrtBox[
+        RowBox[{
+         SuperscriptBox["8", 
+          RowBox[{
+           RowBox[{"-", "2"}], "x"}]], "+", "4"}]], "-", "2"}], ",", "12"}], 
+     "]"}], "//", "ScientificForm"}], ",", 
+   RowBox[{"{", 
+    RowBox[{"x", ",", "1", ",", "10"}], "}"}]}], "]"}]], "Input",
+ CellChangeTimes->{{3.7602543573121824`*^9, 3.7602544343208055`*^9}, {
+  3.7602545735023885`*^9, 3.7602546673307285`*^9}, {3.7602557112771196`*^9, 
+  3.760255713887948*^9}},
+ CellLabel->"In[18]:=",ExpressionUUID->"b5a642b5-7013-42a7-bbe5-da126634b54e"],
+
+Cell[BoxData[
+ RowBox[{"{", 
+  RowBox[{
+   TagBox[
+    InterpretationBox[
+     RowBox[{"\<\"3.90244273517\"\>", "\[Times]", 
+      SuperscriptBox["10", "\<\"-3\"\>"]}],
+     0.00390244273517467060891934471106027601`12.,
+     AutoDelete->True],
+    ScientificForm], ",", 
+   TagBox[
+    InterpretationBox[
+     RowBox[{"\<\"6.10342249558\"\>", "\[Times]", 
+      SuperscriptBox["10", "\<\"-5\"\>"]}],
+     0.00006103422495584600979603589087695134`12.,
+     AutoDelete->True],
+    ScientificForm], ",", 
+   TagBox[
+    InterpretationBox[
+     RowBox[{"\<\"9.53674089033\"\>", "\[Times]", 
+      SuperscriptBox["10", "\<\"-7\"\>"]}],
+     9.5367408903268297692056562946873`12.*^-7,
+     AutoDelete->True],
+    ScientificForm], ",", 
+   TagBox[
+    InterpretationBox[
+     RowBox[{"\<\"1.49011611383\"\>", "\[Times]", 
+      SuperscriptBox["10", "\<\"-8\"\>"]}],
+     1.490116113833650543233247540347`12.*^-8,
+     AutoDelete->True],
+    ScientificForm], ",", 
+   TagBox[
+    InterpretationBox[
+     RowBox[{"\<\"2.32830643640\"\>", "\[Times]", 
+      SuperscriptBox["10", "\<\"-10\"\>"]}],
+     2.3283064364031710175175891639`12.*^-10,
+     AutoDelete->True],
+    ScientificForm], ",", 
+   TagBox[
+    InterpretationBox[
+     RowBox[{"\<\"3.63797880709\"\>", "\[Times]", 
+      SuperscriptBox["10", "\<\"-12\"\>"]}],
+     3.63797880708840422920995016`12.*^-12,
+     AutoDelete->True],
+    ScientificForm], ",", 
+   TagBox[
+    InterpretationBox[
+     RowBox[{"\<\"5.68434188608\"\>", "\[Times]", 
+      SuperscriptBox["10", "\<\"-14\"\>"]}],
+     5.6843418860807207076123`12.*^-14,
+     AutoDelete->True],
+    ScientificForm], ",", 
+   TagBox[
+    InterpretationBox[
+     RowBox[{"\<\"8.88178419700\"\>", "\[Times]", 
+      SuperscriptBox["10", "\<\"-16\"\>"]}],
+     8.8817841970012503512368`12.*^-16,
+     AutoDelete->True],
+    ScientificForm], ",", 
+   TagBox[
+    InterpretationBox[
+     RowBox[{"\<\"1.38777878078\"\>", "\[Times]", 
+      SuperscriptBox["10", "\<\"-17\"\>"]}],
+     1.38777878078144567071471472414543575798506`12.*^-17,
+     AutoDelete->True],
+    ScientificForm], ",", 
+   TagBox[
+    InterpretationBox[
+     RowBox[{"\<\"2.16840434497\"\>", "\[Times]", 
+      SuperscriptBox["10", "\<\"-19\"\>"]}],
+     2.168404344971008867897356166657654672067`12.*^-19,
+     AutoDelete->True],
+    ScientificForm]}], "}"}]], "Output",
+ CellChangeTimes->{{3.7602546101944413`*^9, 3.7602546675718403`*^9}, 
+   3.7602557145470304`*^9},
+ CellLabel->"Out[18]=",ExpressionUUID->"a7dba732-1b71-4587-bf95-054bd7f929a8"]
+}, Open  ]]
+},
+WindowSize->{751, 817},
+WindowMargins->{{Automatic, 186}, {-50, Automatic}},
+Magnification->1.5,
+FrontEndVersion->"11.3 for Microsoft Windows (64-bit) (March 28, 2018)",
+StyleDefinitions->"Default.nb"
+]
+(* End of Notebook Content *)
+
+(* Internal cache information *)
+(*CellTagsOutline
+CellTagsIndex->{}
+*)
+(*CellTagsIndex
+CellTagsIndex->{}
+*)
+(*NotebookFileOutline
+Notebook[{
+Cell[CellGroupData[{
+Cell[580, 22, 478, 10, 131, "Input",ExpressionUUID->"27aa795a-19b1-4daa-99c2-798100ecb41b"],
+Cell[1061, 34, 413, 8, 88, "Output",ExpressionUUID->"3d644f9c-94d5-4c0c-8f97-2c028531f19d"]
+}, Open  ]],
+Cell[CellGroupData[{
+Cell[1511, 47, 200, 3, 43, "Input",ExpressionUUID->"abc47792-61a2-410c-851a-59cf59a6a8bd"],
+Cell[1714, 52, 227, 3, 49, "Output",ExpressionUUID->"6c2dff50-bcc1-49d8-a108-36311f38843c"]
+}, Open  ]],
+Cell[CellGroupData[{
+Cell[1978, 60, 652, 18, 114, "Input",ExpressionUUID->"b5a642b5-7013-42a7-bbe5-da126634b54e"],
+Cell[2633, 80, 2519, 75, 209, "Output",ExpressionUUID->"a7dba732-1b71-4587-bf95-054bd7f929a8"]
+}, Open  ]]
+}
+]
+*)
+
diff --git a/计算方法/labs/2/Lagrange.m b/计算方法/labs/2/Lagrange.m
@@ -0,0 +1,43 @@
+format longE
+syms x
+f = @(x) 1 / (4+x+x.^2);
+N = [4, 8, 16, 128];
+y = (0:500) / 50 - 5;
+polyadd = @(p1, p2) [zeros(1, size(p1,2)-size(p2,2)) p2] + ...
+[zeros(1, size(p2,2)-size(p1,2)) p1];
+
+warning('off', 'all')
+for k = 1:3
+    for kk = 1:2
+        disp(['N = ' num2str(N(k)) ' type' num2str(kk)]);
+        if kk == 1
+            xi = -5 + 10/N(k) * (0:N(k));
+        else
+            xi = -5 * cos((2 * (0:N(k)) + 1)*pi / (2 * N(k) + 2));
+        end
+        p = 0;
+        for i = 0:N(k)
+            l = 1;
+            for j = 0:i-1
+                l = conv(l, [1, -xi(j+1)] / (xi(i+1) - xi(j+1)));
+            end
+            for j = (i+1):N(k)
+                l = conv(l, [1, -xi(j+1)] / (xi(i+1) - xi(j+1)));
+            end
+            p = polyadd(p, f(xi(i+1)) * l);
+        end
+        % p is the result Lagrange polynomial
+        % calculate the max err
+        err = max(abs(arrayfun(f, y) - polyval(p, y)));
+        fprintf("Err = %.12e\n", err);
+        fig = fplot(f, [-5 5]);
+        hold on
+        pval = polyval(p, y);
+        plot(y, pval);
+        title(['N=' num2str(N(k)) ' type=' num2str(kk)]);
+        hold off
+        saveas(fig, ['N' num2str(N(k)) 'type' num2str(kk) '.png']);
+        clf();
+    end
+end
+warning('on', 'all')
diff --git a/计算方法/labs/2/N128type1.png b/计算方法/labs/2/N128type1.png
diff --git a/计算方法/labs/2/N128type2.png b/计算方法/labs/2/N128type2.png
diff --git a/计算方法/labs/2/N16type1.png b/计算方法/labs/2/N16type1.png
diff --git a/计算方法/labs/2/N16type2.png b/计算方法/labs/2/N16type2.png
diff --git a/计算方法/labs/2/N32type1.png b/计算方法/labs/2/N32type1.png
diff --git a/计算方法/labs/2/N32type2.png b/计算方法/labs/2/N32type2.png
diff --git a/计算方法/labs/2/N4type1.png b/计算方法/labs/2/N4type1.png
diff --git a/计算方法/labs/2/N4type2.png b/计算方法/labs/2/N4type2.png
diff --git a/计算方法/labs/2/N64type1.png b/计算方法/labs/2/N64type1.png
diff --git a/计算方法/labs/2/N64type2.png b/计算方法/labs/2/N64type2.png
diff --git a/计算方法/labs/2/N8type1.png b/计算方法/labs/2/N8type1.png
diff --git a/计算方法/labs/2/N8type2.png b/计算方法/labs/2/N8type2.png
diff --git a/计算方法/labs/2/report2.md b/计算方法/labs/2/report2.md
@@ -0,0 +1,73 @@
+# Lab02 Lagrange插值
+
+**古宜民 PB17000002**
+
+**2019.3.10**
+
+### 1. 计算结果
+
+选取插值节点为均匀节点和Chebyshev点对函数 $ f(x)=\frac {1}{4+x+x^2} $ 进行插值，去插值函数和原函数的最大偏差作为近似误差，取不同插值点数N=4，8，16，计算结果如下：
+
+第一组节点：
+N = 4  Err = 6.475332068311e-02
+N = 8  Err = 5.156056628632e-02
+N = 16  Err = 1.071904436126e-01
+
+第二组节点：
+N = 4  Err = 5.437602650073e-02
+N = 8  Err = 1.426092606546e-02
+N = 16  Err = 8.367272096925e-04
+
+#### 函数图像：
+
+其中蓝色为原函数f，橙色为插值函数。
+第一组节点：
+
+![N4type1](.\N4type1.png)
+
+![N8type1](.\N8type1.png)
+
+![N16type1](.\N16type1.png)
+
+第二组节点：
+
+![N4type2](.\N4type2.png)
+
+![N8type2](.\N8type2.png)
+
+![N16type2](.\N16type2.png)
+
+### 2. 程序算法
+
+使用MATLAB。对于给定的N和节点类型，使用两层循环，内层计算插值函数$l_i(x)$，外层计算插值多项式$p(x)=\sum_{i=0}^{n}f(x_i)l_i(x)$。之后在一系列值上计算误差$max\{p(x)-f(x)\}$。
+
+关键代码如下：
+
+```matlab
+        p = 0;
+        for i = 0:N(k)
+            l = 1;
+            for j = 0:i-1
+                l = conv(l, [1, -xi(j+1)] / (xi(i+1) - xi(j+1)));
+            end
+            for j = (i+1):N(k)
+                l = conv(l, [1, -xi(j+1)] / (xi(i+1) - xi(j+1)));
+            end
+            p = polyadd(p, f(xi(i+1)) * l);
+        end
+        % p is the result Lagrange polynomial
+        % calculate the max err
+        err = max(abs(arrayfun(f, y) - polyval(p, y)));
+```
+
+### 3. 结果分析
+
+对于不同次数的插值函数，可见均匀取插值节点时，在N=8时误差最小，N=4、16时误差较大。N=4时误差来自于次数太低，插值函数不能很好地反应原函数性质，偏差较大；而N=16时误差来自于在定义域边缘处出现了Runge现象，误差很大，而中间区域于原函数符合得很好。取Chebyshev节点时，N=4，8，16误差依次减小，N=16时几乎于原函数完全相同。而若取更大的N值，如N=32，则也出现了误差增大的现象。当N更大时，两种插值函数在定义域边缘处都出现了几个数量级的误差。
+
+![N32type2](.\N32type2.png)
+
+相比之下，取Chebyshev节点的插值性质明显好于取均匀节点。N=16时其与原函数几乎完全符合。并且当N过大时产生的偏差也小于均匀节点。
+
+### 4. 小结
+
+由实验结果可见，插值函数的好坏不仅由插值函数的的次数决定，次数过低拟合不好，过高又会在区间端点附近出现较大误差，必须根据作图结果合理选择。另外插值节点的选择也影响结果，Chebyshev节点在端点附近取点较密，中间取点较为稀疏，能得到比均匀取点更稳定的结果、更小的误差。
diff --git a/计算方法/labs/2/report2.pdf b/计算方法/labs/2/report2.pdf
diff --git a/计算方法/labs/3/Simpson.m b/计算方法/labs/3/Simpson.m
@@ -0,0 +1,20 @@
+function [result] = Simpson(f, a, b, N)
+%SIMPSON : Doing Simpson numerical calculus
+%   f: the function to integral
+%   a, b: the intetral interval
+%   N: points number
+h = (b - a) / N;
+result = f(a) + f(b);
+m = N / 2;
+pt1 = 0;
+for i = 0:m-1
+    pt1 = pt1 + f(a + (2 * i + 1) * h);
+end
+result = result + 4 * pt1;
+pt2 = 0;
+for i = 1:m-1
+    pt2 = pt2 + f(a + (2 * i) * h);
+end
+result = result + 2 * pt2;
+result = result * h / 3;
+end