Skip to content

Commit 3add6e4

Browse files
author
petergu
committed
Add Xianjin Chen's computational methods code
1 parent 4aef83c commit 3add6e4

35 files changed

+1289
-0
lines changed

‎计算方法/labs/1/1.cpp

+53
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,53 @@
1+
#include <bits/stdc++.h>
2+
using namespace std;
3+
float f(float x)
4+
{
5+
return sqrt(x * x + 4) - 2;
6+
}
7+
float g(float x)
8+
{
9+
return x * x / (sqrt(x * x + 4) + 2);
10+
}
11+
12+
int main(){
13+
// Problem 1
14+
float x;
15+
x = 1.;
16+
for(int i = 0; i < 10; i++) {
17+
x = x / 8;
18+
cout << setiosflags(ios::scientific) << setprecision(12) << x << "\t" << f(x) << "\t" << g(x) << endl;
19+
}
20+
// Problem 2
21+
double arr[] = {4040.045551380452, -2759471.276702747, -31.64291531266504, 2755462.874010974, 0.0000557052996742893};
22+
double res;
23+
// seq
24+
res = 0.;
25+
for(int i = 0; i < 5; i++) {
26+
res += arr[i];
27+
}
28+
cout << setprecision(7) << "Final " << res << endl;
29+
// rev
30+
res = 0.;
31+
for(int i = 4; i >= 0; i--) {
32+
res += arr[i];
33+
}
34+
cout << "Final " << res << endl;
35+
// posi and nega
36+
res = 0.;
37+
res += arr[0];
38+
res += arr[3];
39+
res += arr[4];
40+
res += arr[1];
41+
res += arr[2];
42+
cout << "Final " << res << endl;
43+
// from small to big
44+
res = 0.;
45+
res += arr[4];
46+
res += arr[2];
47+
res += arr[0];
48+
res += arr[3];
49+
res += arr[1];
50+
cout << "Final " << res << endl;
51+
// Real result is 0.
52+
return 0;
53+
}

‎计算方法/labs/1/1.docx

15.6 KB
Binary file not shown.

‎计算方法/labs/1/cm.nb

+191
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,191 @@
1+
(* Content-type: application/vnd.wolfram.mathematica *)
2+
3+
(*** Wolfram Notebook File ***)
4+
(* http://www.wolfram.com/nb *)
5+
6+
(* CreatedBy='Mathematica 11.3' *)
7+
8+
(*CacheID: 234*)
9+
(* Internal cache information:
10+
NotebookFileLineBreakTest
11+
NotebookFileLineBreakTest
12+
NotebookDataPosition[ 158, 7]
13+
NotebookDataLength[ 6052, 183]
14+
NotebookOptionsPosition[ 5168, 158]
15+
NotebookOutlinePosition[ 5533, 174]
16+
CellTagsIndexPosition[ 5490, 171]
17+
WindowFrame->Normal*)
18+
19+
(* Beginning of Notebook Content *)
20+
Notebook[{
21+
22+
Cell[CellGroupData[{
23+
Cell[BoxData[
24+
RowBox[{
25+
RowBox[{"X", "=",
26+
RowBox[{"{",
27+
RowBox[{"4040.045551380452", ",",
28+
RowBox[{"-", "2759471.276702747"}], ",",
29+
RowBox[{"-", "31.64291531266504"}], ",", "2755462.874010974", ",",
30+
"0.0000557052996742893"}], "}"}]}], "\n"}]], "Input",
31+
CellChangeTimes->{{3.7602515756829224`*^9, 3.760251609278879*^9}, {
32+
3.7602516900422325`*^9, 3.7602516955699854`*^9}},
33+
CellLabel->"In[8]:=",ExpressionUUID->"27aa795a-19b1-4daa-99c2-798100ecb41b"],
34+
35+
Cell[BoxData[
36+
RowBox[{"{",
37+
RowBox[{"4040.045551380452`", ",",
38+
RowBox[{"-", "2.759471276702747`*^6"}], ",",
39+
RowBox[{"-", "31.64291531266504`"}], ",", "2.755462874010974`*^6", ",",
40+
"0.0000557052996742893`"}], "}"}]], "Output",
41+
CellChangeTimes->{
42+
3.760251609675742*^9, {3.7602516904743037`*^9, 3.7602516960183907`*^9}},
43+
CellLabel->"Out[8]=",ExpressionUUID->"3d644f9c-94d5-4c0c-8f97-2c028531f19d"]
44+
}, Open ]],
45+
46+
Cell[CellGroupData[{
47+
48+
Cell[BoxData[
49+
RowBox[{"Total", "[", "X", "]"}]], "Input",
50+
CellChangeTimes->{{3.7602516142985*^9, 3.760251640514464*^9}},
51+
CellLabel->"In[9]:=",ExpressionUUID->"abc47792-61a2-410c-851a-59cf59a6a8bd"],
52+
53+
Cell[BoxData["0.`"], "Output",
54+
CellChangeTimes->{{3.7602516171608987`*^9, 3.760251640820704*^9}, {
55+
3.7602516932206087`*^9, 3.760251698447276*^9}},
56+
CellLabel->"Out[9]=",ExpressionUUID->"6c2dff50-bcc1-49d8-a108-36311f38843c"]
57+
}, Open ]],
58+
59+
Cell[CellGroupData[{
60+
61+
Cell[BoxData[
62+
RowBox[{"Table", "[",
63+
RowBox[{
64+
RowBox[{
65+
RowBox[{"N", "[",
66+
RowBox[{
67+
RowBox[{
68+
SqrtBox[
69+
RowBox[{
70+
SuperscriptBox["8",
71+
RowBox[{
72+
RowBox[{"-", "2"}], "x"}]], "+", "4"}]], "-", "2"}], ",", "12"}],
73+
"]"}], "//", "ScientificForm"}], ",",
74+
RowBox[{"{",
75+
RowBox[{"x", ",", "1", ",", "10"}], "}"}]}], "]"}]], "Input",
76+
CellChangeTimes->{{3.7602543573121824`*^9, 3.7602544343208055`*^9}, {
77+
3.7602545735023885`*^9, 3.7602546673307285`*^9}, {3.7602557112771196`*^9,
78+
3.760255713887948*^9}},
79+
CellLabel->"In[18]:=",ExpressionUUID->"b5a642b5-7013-42a7-bbe5-da126634b54e"],
80+
81+
Cell[BoxData[
82+
RowBox[{"{",
83+
RowBox[{
84+
TagBox[
85+
InterpretationBox[
86+
RowBox[{"\<\"3.90244273517\"\>", "\[Times]",
87+
SuperscriptBox["10", "\<\"-3\"\>"]}],
88+
0.00390244273517467060891934471106027601`12.,
89+
AutoDelete->True],
90+
ScientificForm], ",",
91+
TagBox[
92+
InterpretationBox[
93+
RowBox[{"\<\"6.10342249558\"\>", "\[Times]",
94+
SuperscriptBox["10", "\<\"-5\"\>"]}],
95+
0.00006103422495584600979603589087695134`12.,
96+
AutoDelete->True],
97+
ScientificForm], ",",
98+
TagBox[
99+
InterpretationBox[
100+
RowBox[{"\<\"9.53674089033\"\>", "\[Times]",
101+
SuperscriptBox["10", "\<\"-7\"\>"]}],
102+
9.5367408903268297692056562946873`12.*^-7,
103+
AutoDelete->True],
104+
ScientificForm], ",",
105+
TagBox[
106+
InterpretationBox[
107+
RowBox[{"\<\"1.49011611383\"\>", "\[Times]",
108+
SuperscriptBox["10", "\<\"-8\"\>"]}],
109+
1.490116113833650543233247540347`12.*^-8,
110+
AutoDelete->True],
111+
ScientificForm], ",",
112+
TagBox[
113+
InterpretationBox[
114+
RowBox[{"\<\"2.32830643640\"\>", "\[Times]",
115+
SuperscriptBox["10", "\<\"-10\"\>"]}],
116+
2.3283064364031710175175891639`12.*^-10,
117+
AutoDelete->True],
118+
ScientificForm], ",",
119+
TagBox[
120+
InterpretationBox[
121+
RowBox[{"\<\"3.63797880709\"\>", "\[Times]",
122+
SuperscriptBox["10", "\<\"-12\"\>"]}],
123+
3.63797880708840422920995016`12.*^-12,
124+
AutoDelete->True],
125+
ScientificForm], ",",
126+
TagBox[
127+
InterpretationBox[
128+
RowBox[{"\<\"5.68434188608\"\>", "\[Times]",
129+
SuperscriptBox["10", "\<\"-14\"\>"]}],
130+
5.6843418860807207076123`12.*^-14,
131+
AutoDelete->True],
132+
ScientificForm], ",",
133+
TagBox[
134+
InterpretationBox[
135+
RowBox[{"\<\"8.88178419700\"\>", "\[Times]",
136+
SuperscriptBox["10", "\<\"-16\"\>"]}],
137+
8.8817841970012503512368`12.*^-16,
138+
AutoDelete->True],
139+
ScientificForm], ",",
140+
TagBox[
141+
InterpretationBox[
142+
RowBox[{"\<\"1.38777878078\"\>", "\[Times]",
143+
SuperscriptBox["10", "\<\"-17\"\>"]}],
144+
1.38777878078144567071471472414543575798506`12.*^-17,
145+
AutoDelete->True],
146+
ScientificForm], ",",
147+
TagBox[
148+
InterpretationBox[
149+
RowBox[{"\<\"2.16840434497\"\>", "\[Times]",
150+
SuperscriptBox["10", "\<\"-19\"\>"]}],
151+
2.168404344971008867897356166657654672067`12.*^-19,
152+
AutoDelete->True],
153+
ScientificForm]}], "}"}]], "Output",
154+
CellChangeTimes->{{3.7602546101944413`*^9, 3.7602546675718403`*^9},
155+
3.7602557145470304`*^9},
156+
CellLabel->"Out[18]=",ExpressionUUID->"a7dba732-1b71-4587-bf95-054bd7f929a8"]
157+
}, Open ]]
158+
},
159+
WindowSize->{751, 817},
160+
WindowMargins->{{Automatic, 186}, {-50, Automatic}},
161+
Magnification->1.5,
162+
FrontEndVersion->"11.3 for Microsoft Windows (64-bit) (March 28, 2018)",
163+
StyleDefinitions->"Default.nb"
164+
]
165+
(* End of Notebook Content *)
166+
167+
(* Internal cache information *)
168+
(*CellTagsOutline
169+
CellTagsIndex->{}
170+
*)
171+
(*CellTagsIndex
172+
CellTagsIndex->{}
173+
*)
174+
(*NotebookFileOutline
175+
Notebook[{
176+
Cell[CellGroupData[{
177+
Cell[580, 22, 478, 10, 131, "Input",ExpressionUUID->"27aa795a-19b1-4daa-99c2-798100ecb41b"],
178+
Cell[1061, 34, 413, 8, 88, "Output",ExpressionUUID->"3d644f9c-94d5-4c0c-8f97-2c028531f19d"]
179+
}, Open ]],
180+
Cell[CellGroupData[{
181+
Cell[1511, 47, 200, 3, 43, "Input",ExpressionUUID->"abc47792-61a2-410c-851a-59cf59a6a8bd"],
182+
Cell[1714, 52, 227, 3, 49, "Output",ExpressionUUID->"6c2dff50-bcc1-49d8-a108-36311f38843c"]
183+
}, Open ]],
184+
Cell[CellGroupData[{
185+
Cell[1978, 60, 652, 18, 114, "Input",ExpressionUUID->"b5a642b5-7013-42a7-bbe5-da126634b54e"],
186+
Cell[2633, 80, 2519, 75, 209, "Output",ExpressionUUID->"a7dba732-1b71-4587-bf95-054bd7f929a8"]
187+
}, Open ]]
188+
}
189+
]
190+
*)
191+

‎计算方法/labs/2/Lagrange.m

+43
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,43 @@
1+
format longE
2+
syms x
3+
f = @(x) 1 / (4+x+x.^2);
4+
N = [4, 8, 16, 128];
5+
y = (0:500) / 50 - 5;
6+
polyadd = @(p1, p2) [zeros(1, size(p1,2)-size(p2,2)) p2] + ...
7+
[zeros(1, size(p2,2)-size(p1,2)) p1];
8+
9+
warning('off', 'all')
10+
for k = 1:3
11+
for kk = 1:2
12+
disp(['N = ' num2str(N(k)) ' type' num2str(kk)]);
13+
if kk == 1
14+
xi = -5 + 10/N(k) * (0:N(k));
15+
else
16+
xi = -5 * cos((2 * (0:N(k)) + 1)*pi / (2 * N(k) + 2));
17+
end
18+
p = 0;
19+
for i = 0:N(k)
20+
l = 1;
21+
for j = 0:i-1
22+
l = conv(l, [1, -xi(j+1)] / (xi(i+1) - xi(j+1)));
23+
end
24+
for j = (i+1):N(k)
25+
l = conv(l, [1, -xi(j+1)] / (xi(i+1) - xi(j+1)));
26+
end
27+
p = polyadd(p, f(xi(i+1)) * l);
28+
end
29+
% p is the result Lagrange polynomial
30+
% calculate the max err
31+
err = max(abs(arrayfun(f, y) - polyval(p, y)));
32+
fprintf("Err = %.12e\n", err);
33+
fig = fplot(f, [-5 5]);
34+
hold on
35+
pval = polyval(p, y);
36+
plot(y, pval);
37+
title(['N=' num2str(N(k)) ' type=' num2str(kk)]);
38+
hold off
39+
saveas(fig, ['N' num2str(N(k)) 'type' num2str(kk) '.png']);
40+
clf();
41+
end
42+
end
43+
warning('on', 'all')

‎计算方法/labs/2/N128type1.png

13 KB
Loading

‎计算方法/labs/2/N128type2.png

15.1 KB
Loading

‎计算方法/labs/2/N16type1.png

23.8 KB
Loading

‎计算方法/labs/2/N16type2.png

22.3 KB
Loading

‎计算方法/labs/2/N32type1.png

13.3 KB
Loading

‎计算方法/labs/2/N32type2.png

21.3 KB
Loading

‎计算方法/labs/2/N4type1.png

25.4 KB
Loading

‎计算方法/labs/2/N4type2.png

25.7 KB
Loading

‎计算方法/labs/2/N64type1.png

14.6 KB
Loading

‎计算方法/labs/2/N64type2.png

15.5 KB
Loading

‎计算方法/labs/2/N8type1.png

25 KB
Loading

‎计算方法/labs/2/N8type2.png

25 KB
Loading

‎计算方法/labs/2/report2.md

+73
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,73 @@
1+
# Lab02 Lagrange插值
2+
3+
**古宜民 PB17000002**
4+
5+
**2019.3.10**
6+
7+
### 1. 计算结果
8+
9+
选取插值节点为均匀节点和Chebyshev点对函数 $ f(x)=\frac {1}{4+x+x^2} $ 进行插值，去插值函数和原函数的最大偏差作为近似误差，取不同插值点数N=4，8，16，计算结果如下：
10+
11+
第一组节点：
12+
N = 4 Err = 6.475332068311e-02
13+
N = 8 Err = 5.156056628632e-02
14+
N = 16 Err = 1.071904436126e-01
15+
16+
第二组节点：
17+
N = 4 Err = 5.437602650073e-02
18+
N = 8 Err = 1.426092606546e-02
19+
N = 16 Err = 8.367272096925e-04
20+
21+
#### 函数图像：
22+
23+
其中蓝色为原函数f，橙色为插值函数。
24+
第一组节点：
25+
26+
![N4type1](.\N4type1.png)
27+
28+
![N8type1](.\N8type1.png)
29+
30+
![N16type1](.\N16type1.png)
31+
32+
第二组节点：
33+
34+
![N4type2](.\N4type2.png)
35+
36+
![N8type2](.\N8type2.png)
37+
38+
![N16type2](.\N16type2.png)
39+
40+
### 2. 程序算法
41+
42+
使用MATLAB。对于给定的N和节点类型，使用两层循环，内层计算插值函数$l_i(x)$，外层计算插值多项式$p(x)=\sum_{i=0}^{n}f(x_i)l_i(x)$。之后在一系列值上计算误差$max\{p(x)-f(x)\}$。
43+
44+
关键代码如下：
45+
46+
```matlab
47+
p = 0;
48+
for i = 0:N(k)
49+
l = 1;
50+
for j = 0:i-1
51+
l = conv(l, [1, -xi(j+1)] / (xi(i+1) - xi(j+1)));
52+
end
53+
for j = (i+1):N(k)
54+
l = conv(l, [1, -xi(j+1)] / (xi(i+1) - xi(j+1)));
55+
end
56+
p = polyadd(p, f(xi(i+1)) * l);
57+
end
58+
% p is the result Lagrange polynomial
59+
% calculate the max err
60+
err = max(abs(arrayfun(f, y) - polyval(p, y)));
61+
```
62+
63+
### 3. 结果分析
64+
65+
对于不同次数的插值函数，可见均匀取插值节点时，在N=8时误差最小，N=4、16时误差较大。N=4时误差来自于次数太低，插值函数不能很好地反应原函数性质，偏差较大；而N=16时误差来自于在定义域边缘处出现了Runge现象，误差很大，而中间区域于原函数符合得很好。取Chebyshev节点时，N=4，8，16误差依次减小，N=16时几乎于原函数完全相同。而若取更大的N值，如N=32，则也出现了误差增大的现象。当N更大时，两种插值函数在定义域边缘处都出现了几个数量级的误差。
66+
67+
![N32type2](.\N32type2.png)
68+
69+
相比之下，取Chebyshev节点的插值性质明显好于取均匀节点。N=16时其与原函数几乎完全符合。并且当N过大时产生的偏差也小于均匀节点。
70+
71+
### 4. 小结
72+
73+
由实验结果可见，插值函数的好坏不仅由插值函数的的次数决定，次数过低拟合不好，过高又会在区间端点附近出现较大误差，必须根据作图结果合理选择。另外插值节点的选择也影响结果，Chebyshev节点在端点附近取点较密，中间取点较为稀疏，能得到比均匀取点更稳定的结果、更小的误差。

‎计算方法/labs/2/report2.pdf

292 KB
Binary file not shown.

‎计算方法/labs/3/Simpson.m

+20
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,20 @@
1+
function [result] = Simpson(f, a, b, N)
2+
%SIMPSON : Doing Simpson numerical calculus
3+
% f: the function to integral
4+
% a, b: the intetral interval
5+
% N: points number
6+
h = (b - a) / N;
7+
result = f(a) + f(b);
8+
m = N / 2;
9+
pt1 = 0;
10+
for i = 0:m-1
11+
pt1 = pt1 + f(a + (2 * i + 1) * h);
12+
end
13+
result = result + 4 * pt1;
14+
pt2 = 0;
15+
for i = 1:m-1
16+
pt2 = pt2 + f(a + (2 * i) * h);
17+
end
18+
result = result + 2 * pt2;
19+
result = result * h / 3;
20+
end

0 commit comments

Comments
 (0)