南昌大学数值分析实验报告.doc

- 18 -实 验 报 告实验课程： 数值分析学生名称： 学 号： 专业班级： 2012年 6月1日目 录 一用样条插值法插值的方法生成字体T、5的轮廓3二原子弹爆炸的能量估计14三PageRank算法17 南昌大学实验报告一姓名： 学 号： 专业班级:实验类型：验证 综合设计 创新 实验日期： 实验成绩： 一、 实验目的1、用样条插值的方法生成字体T的轮廓2、C或C+语言用Bezier曲线生成并编写程序二、 实验要求1、 熟悉三次样条插值有关理论，并能将其运用到实际中，加深对理论知识的理解；2、 要求会编程实现Bezier样条曲线，并根据所给数据绘制T的轮廓。三、 主要仪器设备及耗材PC微机Windows 操作系统VS集成开发环境VS0集成开发环境的MSDN四、 实验基本原理和内容- 18 -Times-Roman T的数据X0Y0X1Y1X2Y2X3Y32376202376202371202371202371202373522624143191431914319143014301430143043504350435043504351943519435193532333936339109339109339108339620339620339620507620529602552492552492552492576492576492576492576492570662570662570662570662666266626662666204920492049204922449224492244924860271620183620183620183620237620237620Times-Roman 5的数据X0Y0X1Y1X2Y2X3Y31495971495971495973455973455977615973655993686063686064066953686064066954066953977024066953977023977023826813726763516763516763516763516761426761426763343914267633439334393243832436324343243432428354264442644426744261094201494081494082693723243103242083242083241122643718537185371653714944119661196686906599429942991499087062062024460121012102050282273337833378378123399180399256399256399327381372333422333422288468232491112512112512112512149597149597实验原理：所谓Bezier曲线，是应用于二维应用程序的数学曲线。曲线的定义有四个点：起始点、终止点（也称锚点）以及两个相互分离的中间点。滑动两个中间点，贝塞尔曲线的形状会发生变化。该曲线公式为：假设给出n+1个控制点位置：pk=(xk,yk),0kn，这些坐标点将混合产生下列位置向量P(u)，用来描述p0和pn间逼近Bezier多项式的路径，,0k1，BEZk,n(u)是Bezier混合函数Bernstein多项式：BEZk,n(u)=C(n,k)uk(1-u)n-k，参数C(n,k)是二项式系数：。多数时候，Bezier曲线是一个阶数比控制点少1的多项式,即4点可以生成一个三次曲线。C(n,k)=(n-k+1)/kC(n,k-1),在nk时满足BEZk,n(u)=(1-u)BEZk,n-1(u)+uBEZk-1,n-1(u)，得函数：void binomialCoefficent(GLint n,GLint *c)和void computeBezPnt(GLfloat u,Point2D *bezPnt,GLint nCtrlPnts,Point2D *ctrlPnts,GLint *c)三次Bezier曲线有四个控制点，我们将n=3代入曲线方程中得到BEZ0,3=(1-u)3，BEZ1,3=3u(1-u)2，BEZ2,3=3u2(1-u)，BEZ3,3=u3。得函数void bezier(Point2D *ctrlPnts,GLint nCtrlPnts,GLint nBezCurvePnts)。Bezier曲线生成：1、确定曲线的阶次; 2 计算Bernstein基函数的表达式:生成字体T共分成15段，每段用了三次的代数方程，用Bezier曲线生成当n=3时其中3 把Bezier曲线中的Pk写成分量坐标的形式4 确定一合适的步长；控制t从0到1变化，求出一系列(x,y)坐标点；将其用小线段顺序连接起来。 算法描述：对于二维平面的情况，只有x,y坐标分量，可以给出四点三次Bezier曲线如下的算法描述：输入：阶次，3； 控制顶点：4个，(x0,y0),（x3,y3) begin x=x0 y=y0 moveto (x,y) for t0 to 1 step Dt xB0,3(t)x0B1,3(t)x1B2,3(t)x2B3,3(t)x3 yB0,3(t)y0B1,3(t)y1B2,3(t)y2B3,3(t)y3 lineto (x,y) endfor end三次Bezier曲线例子：对于平面上给定的7个控制点坐标分别为:A（100,300）,B（120,200）,C（220,200）,D（270,100）,E（370,100）, F（420,200）,G（420,300）。画出其三次Bezier曲线。通过以上的准备我们开始编写主函数，再输入T的实验数据，生成图像。五、实验程序及运行结果 实验程序代码如下/ sy1.cpp : Defines the entry point for the application.#include stdafx.h#include resource.h#define MAX_LOADSTRING 100/ Global Variables:HINSTANCE hInst;/ current instanceTCHAR szTitleMAX_LOADSTRING;TCHAR szWindowClassMAX_LOADSTRING;/ The title bar textATOMMyRegisterClass(HINSTANCE hInstance);BOOLInitInstance(HINSTANCE, int);LRESULT CALLBACKWndProc(HWND, UINT, WPARAM, LPARAM);LRESULT CALLBACKAbout(HWND, UINT, WPARAM, LPARAM);int APIENTRY WinMain(HINSTANCE hInstance, HINSTANCE hPrevInstance, LPSTR lpCmdLine, int nCmdShow) MSG msg;HACCEL hAccelTable; LoadString(hInstance, IDS_APP_TITLE, szTitle, MAX_LOADSTRING);LoadString(hInstance, IDC_SY1, szWindowClass, MAX_LOADSTRING);MyRegisterClass(hInstance); if (!InitInstance (hInstance, nCmdShow) return FALSE;hAccelTable = LoadAccelerators(hInstance, (LPCTSTR)IDC_SY1);while (GetMessage(&msg, NULL, 0, 0) if (!TranslateAccelerator(msg.hwnd, hAccelTable, &msg) TranslateMessage(&msg);DispatchMessage(&msg); return msg.wParam;ATOM MyRegisterClass(HINSTANCE hInstance)WNDCLASSEX wcex;wcex.cbSize = sizeof(WNDCLASSEX); wcex.style= CS_HREDRAW | CS_VREDRAW;wcex.lpfnWndProc= (WNDPROC)WndProc;wcex.cbClsExtra= 0;wcex.cbWndExtra= 0;wcex.hInstance= hInstance;wcex.hIcon= LoadIcon(hInstance, (LPCTSTR)IDI_SY1);wcex.hCursor= LoadCursor(NULL, IDC_ARROW);wcex.hbrBackground= (HBRUSH)(COLOR_WINDOW+1);wcex.lpszMenuName= (LPCSTR)IDC_SY1;wcex.lpszClassName= szWindowClass;wcex.hIconSm= LoadIcon(wcex.hInstance, (LPCTSTR)IDI_SMALL); return RegisterClassEx(&wcex);BOOL InitInstance(HINSTANCE hInstance, int nCmdShow) HWND hWnd; hInst = hInstance; / Store instance handle in our global variable hWnd = CreateWindow(szWindowClass, szTitle, WS_OVERLAPPEDWINDOW, CW_USEDEFAULT, 0, CW_USEDEFAULT, 0, NULL, NULL, hInstance, NULL); if (!hWnd) return FALSE; ShowWindow(hWnd, nCmdShow); UpdateWindow(hWnd); return TRUE;LRESULT CALLBACK WndProc(HWND hWnd, UINT message, WPARAM wParam, LPARAM lParam)int wmId, wmEvent,i;PAINTSTRUCT ps;HDC hdc;TCHAR szHelloMAX_LOADSTRING;LoadString(hInst, IDS_HELLO, szHello, MAX_LOADSTRING); float p(int point0,int point1,int point2,int point3 ,float t);switch (message) case WM_COMMAND:wmId = LOWORD(wParam); wmEvent = HIWORD(wParam); switch (wmId)case IDM_ABOUT:DialogBox(hInst, (LPCTSTR)IDD_ABOUTBOX, hWnd, (DLGPROC)About); break;case IDM_EXIT: DestroyWindow(hWnd); break;default: return DefWindowProc(hWnd, message, wParam, lParam);break;case WM_PAINT:hdc = BeginPaint(hWnd, &ps);/ TODO: Add any drawing code here.RECT rt;GetClientRect(hWnd, &rt);/DrawText(hdc, szHello, strlen(szHello), &rt, DT_CENTER); for(i=0;i1000;i+)SetPixel(hdc,800-Pt(237,237,237,237,0.001*i),680-Pt(620,620,120,120,0.001*i),RGB(250,0,0); SetPixel(hdc,800-Pt(237,237,226,143,0.001*i),680-Pt(120,35 ,24 ,19 ,0.001*i),RGB(250,0,0); SetPixel(hdc,800-Pt(143,143,143,143,0.001*i),680-Pt(19 ,19,0,0,0.001*i),RGB(250,0,0);SetPixel(hdc,800-Pt(143,143,435,435,0.001*i),680-Pt( 0, 0, 0, 0,0.001*i),RGB(250,0,0);SetPixel(hdc,800-Pt(435,435,435,435,0.001*i),680-Pt( 0, 0,19,19,0.001*i),RGB(250,0,0); SetPixel(hdc,800-Pt(435,353,339,339,0.001*i),680-Pt(19,23,36,109,0.001*i),RGB(250,0,0); SetPixel(hdc,800-Pt(339,339,339,339,0.001*i),680-Pt(109,108,620,620,0.001*i),RGB(250,0,0); SetPixel(hdc,800-Pt(339,507,529,552,0.001*i),680-Pt(620,620,602,492,0.001*i),RGB(250,0,0); SetPixel(hdc,800-Pt(552,552,576,576,0.001*i),680-Pt(492,492,492,492,0.001*i),RGB(250,0,0); SetPixel(hdc,800-Pt(576,576,570,570,0.001*i),680-Pt(492,492,662,662,0.001*i),RGB(250,0,0); SetPixel(hdc,800-Pt(570,570,6,6,0.001*i),680-Pt(662,662,662,662,0.001*i),RGB(250,0,0); SetPixel(hdc,800-Pt(6,6,0,0,0.001*i),680-Pt(662,662,492,492,0.001*i),RGB(250,0,0); SetPixel(hdc,800-Pt(0,0,24,24,0.001*i),680-Pt(492,492,492,492,0.001*i),RGB(250,0,0); SetPixel(hdc,800-Pt(24,48,71,183,0.001*i),680-Pt(492,602,620,620,0.001*i),RGB(250,0,0); SetPixel(hdc,800-Pt(183,183,237,237,0.001*i),680-Pt(620,620,620,620,0.001*i),RGB(250,0,0);EndPaint(hWnd, &ps);break;case WM_DESTROY:PostQuitMessage(0);break;default:return DefWindowProc(hWnd, message, wParam, lParam); return 0;/ Mesage handler for about box.LRESULT CALLBACK About(HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam)switch (message)case WM_INITDIALOG:return TRUE; case WM_COMMAND:if (LOWORD(wParam) = IDOK | LOWORD(wParam) = IDCANCEL) EndDialog(hDlg, LOWORD(wParam);return TRUE;break; return FALSE;float p(int point0,int point1,int point2,int point3 ,float t) float pp;pp= point0* (1 - t) * (1 - t) * (1 - t) + point1* 3 * t * (1 - t) * (1 - t) + point2 * 3 * (t * t) * (1 - t) + point3 * (t * t * t);return pp;若要得到5字，只需改for循环那部分：for(i=0;i1000;i+)SetPixel(hdc,500-Pt(-149,-149,-149,-345,0.001*i),700-Pt(597,597,597,597,0.001*i),RGB(0,255,0);SetPixel(hdc,500-Pt(-345,-365,-365,-368,0.001*i),700-Pt(597,597,599,606,0.001*i),RGB(0,255,0);SetPixel(hdc,500-Pt(-368,-406,-368,-406,0.001*i),700-Pt(606,695,606,695,0.001*i),RGB(0,255,0);SetPixel(hdc,500-Pt(-406,-397,-406,-397,0.001*i),700-Pt(695,702,695,702,0.001*i),RGB(0,255,0);SetPixel(hdc,500-Pt(-397,-382,-372,-351,0.001*i),700-Pt(702,681,676,676,0.001*i),RGB(0,255,0);SetPixel(hdc,500-Pt(-351,-351,-351,-142,0.001*i),700-Pt(676,676,676,676,0.001*i),RGB(0,255,0);SetPixel(hdc,500-Pt(-142,-33,-142,-33,0.001*i),700-Pt(676,439,676,439,0.001*i),RGB(0,255,0);SetPixel(hdc,500-Pt(-33,-32,-32,-32,0.001*i),700-Pt(439,438,436,434,0.001*i),RGB(0,255,0);SetPixel(hdc,500-Pt(-32,-32,-35,-44,0.001*i),700-Pt(434,428,426,426,0.001*i),RGB(0,255,0);SetPixel(hdc,500-Pt(-44,-74,-109,-149,0.001*i),700-Pt(426,426,420,408,0.001*i),RGB(0,255,0);SetPixel(hdc,500-Pt(-149,-269,-324,-324,0.001*i),700-Pt(408,372,310,208,0.001*i),RGB(0,255,0);SetPixel(hdc,500-Pt(-324,-324,-264,-185,0.001*i),700-Pt(208,112,37,37,0.001*i),RGB(0,255,0);SetPixel(hdc,500-Pt(-185,-165,-149,-119,0.001*i),700-Pt(37,37,44,44,0.001*i),RGB(0,255,0);SetPixel(hdc,500-Pt(-119,-86,-65,-42,0.001*i),700-Pt(66,90,99,99,0.001*i),RGB(0,255,0);SetPixel(hdc,500-Pt(-42,-14,-0,-0,0.001*i),700-Pt(99,99,87,62,0.001*i),RGB(0,255,0);SetPixel(hdc,500-Pt(-0,-0,-46,-121,0.001*i),700-Pt(62,24,0,0,0.001*i),RGB(0,255,0);SetPixel(hdc,500-Pt(-121,-205,-282,-333,0.001*i),700-Pt(0,0,27,78,0.001*i),RGB(0,255,0);SetPixel(hdc,500-Pt(-333,-378,-399,-399,0.001*i),700-Pt(78,123,180,256,0.001*i),RGB(0,255,0);SetPixel(hdc,500-Pt(-399,-399,-381,-333,0.001*i),700-Pt(256,327,372,422,0.001*i),RGB(0,255,0);SetPixel(hdc,500-Pt(-333,-288,-232,-112,0.001*i),700-Pt(422,468,491,512,0.001*i),RGB(0,255,0);SetPixel(hdc,500-Pt(-112,-112,-149,-149,0.001*i),700-Pt(512,512,597,597,0.001*i),RGB(0,255,0);运行结果显示为：五、 实验总结和体会 图形类实验，第一次，也是第一次成功。希望能在以后的编程生涯中打下坚实的基础。 南昌大学实验报告二姓名： 学 号： 专业班级: 实验类型：验证 综合设计 创新 实验日期： 实验成绩： 一、实验目的用C语言或C+语言编程实现原子弹爆炸的能量估计。二、实验基本原理和内容1945年7月16 日，美国科学家在新墨西哥州Los Alamos 沙漠试爆了世界上第一颗原子弹，这一事件令全球震惊。但在当时有关原子弹爆炸的任何资料都是保密的，而很多其他国家的科学家非常想知道这次爆炸的威力有多大。两年之后，美国政府首次公开了这次爆炸的录像带，而其他的数据和资料仍然不被外界所知，英国物理学家G.I.Taylor(1886-1975)通过研究原子弹爆炸的录像带，建立数学模型对爆炸所释放的能量进行了估计，得到估计值与若干年后正式公布的爆炸能量21kt相当接近（1kt为1千吨TNT炸药的爆炸能量）。Taylor是如何根据爆炸录像估计的呢？主要是通过测量爆炸形成的“蘑菇云”半径来进行估计的。因为爆炸产生的冲击波从中心点向外传播，爆炸的能量越大，在相同时间内冲击波传播得越远，蘑菇云的半径就越大，Taylor通过研究录像带，测量了从爆炸开始的不同时刻t所对应的蘑菇云半径r(t)，如下表所示tr(t)tr(t)tr(t)tr(t)tr(t)0.1011.10.8034.21.5044.43.5361.115.0106.50.2419.90.9436.31.6546.03.8062.925.0130.00.3825.41.0838.91.7946.94.0764.334.0145.00.5228.81.2241.01.9348.74.3465.653.0175.00.6631.91.3642.83.2659.04.6197.362.0185.0*t的单位为ms，r的单位为m。然后通过量纲分析法建立了蘑菇云半径r与时间t和爆炸能量E的关系，利用上述数据最后求出了爆炸的能量。三、主要仪器设备及耗材PC机一台PC微机Windows 操作系统Microsoft Visual Studio 6.0集成开发环境一、 实验程序及运行结果 实验程序代码如下：#include stdio.h#include math.h#define e 2.718281828459 int i,m=24;double neiji(double *g,double *h) double sum=0;for(i=0;i=m;i+)sum+=gi*hi;return sum;main()double t=0;double A,a,b,p=1.25;double a0,a1,a2,b0,b1,c,E;double s025;double s125=0.1,0.24,0.38,0.52,0.66,0.8,0.94,1.08,1.22,1.36,1.5,1.65,1.79,1.93,3.26,3.53,3.8,4.07,4.34,4.61,15,25,34,53,62;double f25=11.1,19.9,25.4,28.8,31.9,34.2,36.3,38.9,41,42.8,44.4,46,46.9,48.7,59,61.1,62.9,64.3,65.6,97.3,106.5,130,145,175,185;double k25;for(i=0;i=24;i+)s0i=1;s1i=0.001*s1i; /单位由毫秒变为秒 a0=neiji(s0,s0);a1=neiji(s0,s1);a2=neiji(s1,s1);b0=neiji(s0,f);b1=neiji(s1,f);printf(s0,s0)=%lfn,neiji(s0,s0);printf(s0,s1)=(s1,s0)=%lfn,neiji(s0,s1);printf(s1,s1)=%lfn,neiji(s1,s1);printf(s0,f)=%lfn,neiji(s0,f);printf(s1,f)=%lfn,neiji(s1,f);b=(a0*b1-a1*b0)/(a2*a0-a1*a1);A=(b0-a2*b)/a0; printf(b=%lfn,b);a=pow(e,A);printf(A=%lf,b=%lfn,A,log(b); for(i=0;i=25;i+)ki=5*log(fi)-2*log(s1i); c=neiji(k,s0)/a0;E=p*pow(e,c);printf(E=%e,爆炸能量=%ektn,E,E/(4.184*pow(10,12);return 0;结果及截图：5、 实验总结及体会 以前只觉得数值分析这门课程只是一门十分繁琐的课程，高数线代等代数内容都攘括在内，而且都是一些非常复杂的代数算法公式，学习起来十分不情愿。但通过此次试验让我了解到，尽管课程十分繁琐，但同样能解决非常繁琐的问题，像本实验中原子弹爆炸的能量估算问题，受益匪浅。 南昌大学实验报告三姓名： 学 号： 专业班级: 实验类型：验证 综合设计 创新 实验日期： 实验成绩： 一、 实验目的用C语言或C+编程实现Google 的PageRank算法二、 实验基本原理和内容PageRank技术的基本原理是：如果网页A链接到网页B，则认为“网页A投了网页B一票”，而且如果网页A是级别高的网页，则网页B的级别也相应地高。要计算PageRank，可假设一个随机上网“冲浪”的过程，即每次看完当前网页后，有两种选择：（1） 在当前网页中随机选一个超链接进入下一个网页；（2） 随机地新开一个网页。这在数学上称为马尔可夫过程，若这样的随机“冲浪”一直进行下去，某个网页被访问到的极限概率就是它的PageRank。设p为选择当前网页上链接的概率（比如，p=0.85），则1-p为不选当前网页的链接而随机打开一个网页的概率。若当前网页是网页j，则如何计算下一步浏览到达网页i的概率（网页j到i的转移概率）？它有两种可能性：（1） 若网页i在网页j的链接上，其概率为p1/cj+(1-p) 1/n；（2）若网页i不在网页j的链接上，其概率为(1-p) 1/n由于网页i是否在网页j的链接上由gij决定，网页j到i的转移概率为： 应注意到的是,若cj=0意味着gij=0，上式改为aij=1/n。任意两个网页之间的转移概率形成了一个转移矩阵A=（aij），设矩阵D为各个网页出度的倒数（若没有出度，设为1）构成的n阶对角阵，e为全是1的n维向量，则:设表示某时刻k浏览网页i的概率，向量x(k)表示当前时刻浏览各网页的概率分布。那么下一时刻浏览到网页i的概率为，此时浏览各网页的概率分布为x(k+1)=Ax(k)当这个过程无限进行下去，达到极限情况，即网页访问概率x(k)收敛到一个极限值，这个极限向量x为各网页的PageRank，它满足Ax=x,且三、 主要仪器设备及耗材PC机一台四、 实验程序及运行结果PageRank.h#include #include #include #include using namespace std; / use graph store webpage, weight representlink timesclass Node public: explicit Node(string name, double pr = 1):name_(name), page_rank_(pr) Node() linkin_nodes_.clear(); void InsertLinkdInNode(Node* node) /如果没有链接 if (linkin_nodes_.find(node) = linkin_nodes_.end() linkin_nodes_.insert(node); node-InsertLinkOutNode(this); void InsertLinkOutNode(Node* node) /如果没有链接 if (linkout_nodes_.find(node) = linkout_nodes_.end() linkout_nodes_.insert(node); double GetPageRank() return page_rank_; void SetPageRank(double pr) page_rank_ = pr; double CalcRank() double pr = 0; set:const_iterator citr = linkin_nodes_.begin(); for (; citr != linkin_nodes_.end(); +citr) Node * node = *citr; pr += node-GetPageRank()/node-GetOutBoundNum(); return pr; size_t GetOutBoundNum() return linkout_nodes_.size(); size_t GetInBoundNum() return linkin_nodes_.size(); void PrintNode() cout Node: name_ s pagerank is: page_rank_ endl; private: string name_; set linkin_nodes_; set linkout_nodes_; double page_rank_; class PageRank public: PageRank(double q=0.85); PageRank(void); void Calc(vector & nodes, int n); double Calc(Node* node); void PrintPageRank(v