c++ 基础向量类【数值计算01】

 这学期的《数值计算》课程挺有意思的，将数学公式转成代码是一个工程师必不可少的能力。

而《数值计算》课就是完成这么一步的课程，实现的语言是c++，主要是自己实现向量，矩阵这些基础类，然后实现一些解方程的算法。《数值计算》每周3个学时，每周都会有代码任务。所以每周在做的时候写点笔记发不出来，给自己更多思考的机会。

作业要求：

第一部分
Vector.h和Vector.cpp两个文件里面写了一个关于向量的模板类，
请大家查阅C++类和向量计算的有关知识，并完成如下任务：
1. 建立一个C++ 空项目，并加入这两个文件。（1分）
2. C++项目里面创建一个main.cpp文件，并在里面写main函数，
main函数里面自行定义两个大小为100的向量，
向量的元素值随机(要用到生成随机数的函数，如何使用自行百度)。（1分）
3. 利用Vector类里面的运算符重载函数，在main函数里面实现这两个向量的加、减、乘操作。（1分）
4. 在Vector类里面编写向量叉乘和计算向量2范数的函数，
并在main函数里面定义两个大小为3的向量，计算两个向量的叉乘和2范数。（1分）
5. 所有的结果均要能够输出在txt文档里面。（1分）

第二部分
用Matlab实现圆周率的计算

先把代码链接奉上

我的实现过程：

首先声明下vector 模板类的结构

1. template <typename T> 
2. class Vector{
3. private:
4.     int num;
5.     T *data;
6. }

1、实现随机数

c语言和c++有不同的实现随机数的方法。

c语言是在 stdlib.h 里面的 rand()函数

c++是 在 #include <random>里面的 std::random_device 类，该类重载了函数符号（）

使用demo如下,定义的对象就想是函数那样使用。

1. std::random_device rd;
2. rd();

所以在c++里面就有两种选择了。

原先想到用rand函数

1. srand(seed);
2.         for (int i = 0; i < n; i++) {
3.             data[i] = T(rand());
4.         }    

但是问题来了 seed怎么设置，seed设置成一样的话，会有两个问题，第一，不同对象的随机数都一样，程序上下文定义的对象 数值都一样，第二，下一次再运行时，对象数值也一样。一开始想用当前时间作为随机种子，可以解决第一个问题。但是第二个问题又有上下文的时间的差值在小数点后很多位数，而随机种子需要的是整数，在类型转换过程中，把两个数的差异丢失了。后来想采用c++的random_device 来产生随机数，能够完美的解决上述来那个个问题，但是random头文件里面包含了标准库的vector容器，跟我们的自己的类名产生了冲突，不得已放弃。

最后解决，既然要求随机种子在不同对象定义的时候不一样，那么就来个不一样的整数就行了，于是想到了地址，每隔对象的成员变量的地址肯定是不一样的，虽然在上下文定义的两个对象，他们的成员变量的地址很接近，但地址是四个字节的，整形在64位机子上也是4 byte的，所以差异还是很够分辨出来的，我把来那个个对象的num变量的地址打印出来看看：

        seed 7338088 & num 006FF868 size 4
        seed 7338072 & num 006FF858 size 4

总而言之，有差别就行了。。好了，随机数讲完了。

2、运算符的重载

运算符的重载就比较多了，挑一个最重要的来讲，[] 中括号运算符。

T & operator [] (int i) { return data[i]; }

重载后我们可以直接通过下标来获取元素的值，但是获取到了元素的值，比如 int a = vector_1[2] 就能获取到第三个数据的值。特别注意前 T & operator 这个& 引用符号，之后加上引用符号，才能够直接对数据进行赋值。比如 vector_1[2] = 10。

3、向量的叉乘

这里我们只讨论长度为3的向量的叉乘。在线性代数中向量的叉乘是运算如下：

代码实现如下：

1. template <typename T>
2. void cross_product(Vector<T> &vec1, Vector<T> &vec2, Vector<T> &vec){
3.     if (vec1.size() != 3 || vec2.size() != 3 || vec.size() != 3) {
4.         cerr << "input vector size != 3" << endl;
5.     }
6.         vec[0] = vec1[1] * vec2[2] - vec1[2] * vec2[1];
7.         vec[1] = vec1[2] * vec2[0] - vec1[0] * vec2[2];
8.         vec[2] = vec1[0] * vec2[1] - vec1[1] * vec2[0];
9. }

4、二范数

  二范数就是向量各个元素之和再开根号，需要注意的是conj()是取共轭的，谁啊让我们这里的元素并不是复数，老师让加上就加上吧

1. template <typename T>
2. double norm2(Vector<T> &vec) {
3.     T sum2 = 0.0;
4.     for (int i = 0; i < vec.size(); i++) {
5.         sum2 += conj(vec[i]) * vec[i];
6.     }
7.     
8.     return sqrt(sum2);
9. }

5、蒙特卡洛方法求pi

思路很简单，定义一个圆和一个正方形，圆的直径等于正方形的边长，圆的面积是 s = pi *r*r ，正方形的面积是 s =2 r * 2r。所以二者一比就得出了pi / 4。我们往正方形里面随即撒点，用点数来代表面积，统计下落在圆内的点数，用这个数字除以所有点数，就能近似的求出面积比，如果点数越多，比例就越接近于真实的比例。代码实现中有几个关键点，一个点的横坐标x和纵坐标y服从在（-r，+r）的均匀分布。第二个是，怎么统计点在圆内，思路也很简单，根据圆的定义，计算点的坐标和原点坐标的距离（假设我们将原点设为圆心），如果距离小于半径，那么统计数就加一，对每一个点都计算一遍，落在圆内的点数就出来了。

 先上结果：

正态分布  10000 点 统计出的结果是 3.1328

代码如下：

1. import numpy as np
2. import matplotlib.pyplot as plt
3. from matplotlib.patches import Circle
4. # 投点次数
5. n = 10000
6.  
7. # 圆的信息
8. r = 1.0         # 半径
9. a, b = (0., 0.) # 圆心
10. # 正方形区域边界
11. x_min, x_max = a-r, a+r
12. y_min, y_max = b-r, b+r
13. # 在正方形区域内随机投点
14. x = np.random.uniform(x_min, x_max, n) # 均匀分布
15. y = np.random.uniform(y_min, y_max, n)
16.  
17. # 计算 点到圆心的距离
18. d = np.sqrt((x-a)**2 + (y-b)**2)
19. # 统计 落在圆内的点的数目
20. res = sum(np.where(d < r, 1, 0))
21. # 计算 pi 的近似值（Monte Carlo方法的精髓：用统计值去近似真实值）
22. pi = 4 * res / n
23. print('pi: ', pi)
24. # 画个图看看
25. fig = plt.figure() 
26. axes = fig.add_subplot(111)
27. axes.plot(x, y,'yo',markersize = 0.5, label='The 333333333333333333333red data')
28. plt.axis('equal') # 防止图像变形
29.  
30. circle = Circle(xy=(a,b), radius=r,alpha=0.5)
31. axes.add_patch(circle)
32. plt.savefig("uniform_10000.jpg")
33. plt.show()

over