数值计算实习.
北航数值分析计算实习题目二 矩阵QR分解
数值分析实习二院（系）名称航空科学与工程学院专业名称动力工程及工程热物理学号SY0905303学生姓名解立垚1. 题目试用带双步位移QR 的分解法求矩阵A=[a ij ]10*10的全部特征值，并对其中的每一个实特征值求相应的特征向量。
已知()sin 0.50.2,1.5cos 1.2,ij i j i j a i j i j ⎧⎫+≠⎪⎪=⎨⎬+=⎪⎪⎩⎭(),1,2,...,10i j =。
说明：1、求矩阵特征值时，要求迭代的精度水平为1210ε-=。
2、打印以下内容：算法的设计方案；全部源程序(要求注明主程序和每个子程序的功能)； 矩阵A 经过拟上三角话之后所得的矩阵()1n A -；对矩阵()1n A-进行QR 分解方法结束后所得的矩阵；矩阵A 的全部特征值()(),1,2,......10i i iR I i λ=，和A 的相应于实特征值的特征向量；其中()(),.i e i m i R R I I λλ==如果i λ是实数，则令0.i I =3、采用e 型输出数据，并且至少显示12位有效数字。
2. 算法设计方案本题采用带双步位移的QR 分解方法。
为了使程序简洁，自定义类Xmatrix ，其中封装了所需要的函数方法。
在Xmatrix 类中封装了运算符重载的函数，即定义了矩阵的加、减、乘、除、数乘运算及转置运算(T())。
同时为了避免传递数组带来的额外内存开销，使用引用(&)代替值传递，以节省内存空间，避免溢出.（1）此程序的主要部分为Xmatrix 中的doubleQR()方法，具体如下：Step1：使用矩阵拟上三角化的算法将A 化为拟上三角阵A （n-1）(此处调用Xmatrix 中的preQR()方法)Step2：令121,,10k m n ε-===， 其中k 为迭代次数。
Step3：如果,1m m a ε-≤，则得到A 的一个特征值,m m a ，令1m m =-，goto Step4;否则goto Step5.Step4： 如果1m =，则得到A 的一个特征值11a ，goto Step11;如果0m =，则goto Step11;如果1m >，则goto Step3;Step5(Step6)：如果2m =，则得到A 的两个特征值12s s 和（12s s 和为右下角两阶子阵对应的特征方程21,1,()det 0m m m m a a D λλ---++=的两个根。
《数值计算方法实习》教学大纲
《数值计算方法实习》教学大纲Numerical Computation Method Practice适用本科四年制信息与计算科学专业（2周 2学分）一、课程的目的和任务本课程的授课对象是信息与计算科学专业本科生，属信息与计算科学专业公共基础课。
数值计算方法是一门专门研究各种数学问题近似解法的课程，它是一门与计算机应用密切结合的实用性很强的数学课程。
在数值计算方法课程中，讲授了各种数学问题的近似解法，这些近似解法的计算量很大，只有利用计算机计算，这些解法才具有实用意义。
因而上机实习，掌握这些近似解法的计算机实现是数值计算方法课程学习的一个重要环节。
本课程实习的主要目的是通过科学计算语言MA TLAB的学习，利用MA TLAB求解各种数学问题的近似解，使学生对数值计算方法课程所学的各种近似解法能在计算机上实现，提高学生对数值计算方法课程讲授的各种数学问题近似解法的理解和掌握。
通过本实践环节，要求学生初步掌握MATLAB的使用方法，掌握利用MATLAB求解各种数学问题近似解的算法，通过上机实践，提高学生对各种数学问题近似解法的实际运用能力，并能应用所学的方法解决一些较简单的实际问题。
二、课程的基本要求和特点本课程是一门既有系统理论又有较强实践性的技术基础课，学习本课程需坚持理论联系实际的学风，必须在学习数值计算方法课程讲授的各种数学问题近似解法的基础上，动手编写一些简单的MA TLAB程序，利用MATLAB来实现求解各种数学问题的近似解；同时要注意数学软件的使用原理及使用方法。
本课程是一门实用性很强的应用数学课程。
三、本课程与其它课程的联系本课程实习是对前期《数值计算方法》课程的巩固，数值计算方法课程涉及面较宽，必须先修课程为《数学分析》、《高等代数》、《常微分方程》、《计算机应用基础》、《数值计算方法》。
四、课程的主要内容1 数学软件MATLAB教学要求：了解：MA TLAB的基本特点，MATLAB的启动方法和工作界面，MATLAB数值计算，MATLAB程序设计，MATLAB绘图。
数值计算方法实习报告
实习报告实习单位：XX大学计算中心实习时间：2023年1月1日至2023年1月31日实习内容：数值计算方法一、实习背景及目的随着科技的不断发展，数值计算方法在工程、物理、化学、生物学等领域发挥着越来越重要的作用。
为了更好地将所学知识应用于实际问题，提高自己的实践能力，我选择了数值计算方法作为实习内容。
本次实习的主要目的是：1. 加深对数值计算方法的理解，掌握基本的数值计算方法及其应用。
2. 提高编程能力，熟练运用C语言进行数值计算程序的设计与实现。
3. 学会分析并解决实际问题，将所学知识运用到实际项目中。
二、实习过程及收获1. 实习前期，我首先学习了数值计算方法的基本理论，包括误差分析、插值法、数值积分、常微分方程数值解等。
通过理论的学习，我对数值计算方法有了更深入的了解。
2. 在实习过程中，我使用C语言编写了一系列数值计算程序，包括求解方程的迭代法、高斯消去法、牛顿法等。
这些程序可以帮助我更好地理解数值计算方法的理论，并提高我的编程能力。
3. 针对实际问题，我运用所学知识进行了解决。
例如，我使用数值积分方法计算了函数在一个区间上的定积分，使用常微分方程数值解方法求解了一个实际物理问题。
这些实践经历使我更加熟悉了数值计算方法在实际问题中的应用。
4. 实习期间，我还参加了计算中心组织的讲座和讨论，与其他实习生交流心得，共同解决问题。
这使我受益匪浅，不仅提高了自己的实际操作能力，还拓宽了知识面。
三、实习总结通过本次实习，我对数值计算方法有了更全面的认识，掌握了基本的数值计算方法及其编程实现。
同时，我的编程能力和解决实际问题的能力也得到了很大提高。
此外，我还学会了如何将所学知识应用于实际项目，为将来的工作打下了坚实基础。
在今后的工作中，我将继续努力学习数值计算方法及相关知识，不断提高自己的实践能力。
同时，我也将把所学知识运用到实际工作中，为公司的发展做出贡献。
最后，感谢计算中心给我提供了一次宝贵的实习机会，使我受益匪浅。
数值方法计算实习题
数值⽅法计算实习题数值⽅法计算实习题⼀、下表给出了飞⾏中鸭⼦的上部形状的节点数据，试⽤三次样条插值函数（⾃然边界条件）和20次Lagrange 插值多项式对数据进⾏插值。
⽤图⽰出给定的数据，以及()s x 和20()L x 。
12 12.6 13.0 13.3];>> y=[1.3 1.5 1.85 2.1 2.6 2.7 2.4 2.15 2.05 2.1 2.25 2.3 2.25 1.95 1.4 0.9 0.7 0.6 0.5 0.4 0.25]; %（1）三次样条插值法xi=0.9:0.01:13.3;yi=interp1(x,y,xi,'spline'); >> xi=0.9:0.01:13.3;yi=interp1(x,y,xi,'spline'); >> title('试验⼀--三次样条插值图⽰')024********试验⼀--三次样条插值图⽰>> pp=spline(x,y)pp =form: 'pp'breaks: [1x21 double]coefs: [20x4 double]pieces: 20order: 4dim: 1>> pp.coefsans =0.7735 -0.9995 0.7760 1.3000 0.7735 -0.0714 0.3477 1.5000 -2.7894 1.3209 1.0974 1.8500 -0.4585 -0.3528 1.2910 2.10000.4489 -1.0405 0.5944 2.6000 0.1738 -0.5018 -0.0225 2.70000.0783 -0.0325 -0.5033 2.40001.3141 0.0850 -0.47712.1500 -1.5812 1.2676 -0.0713 2.0500 0.0431 -0.1555 0.2623 2.1000 -0.0047 -0.0261 0.0808 2.2500 -0.0245 -0.0401 0.0146 2.3000 0.0175 -0.1135 -0.1390 2.2500 -0.0128 -0.0505 -0.3358 1.9500 -0.0201 -0.1003 -0.5319 1.4000 1.2094 -0.1485 -0.7310 0.9000 -0.8279 0.9400 -0.4935 0.7000 0.0122 -0.0535 -0.1389 0.6000 -0.2960 -0.0316 -0.1900 0.5000 -0.2960 -0.3867 -0.3573 0.4000 所以所得⽅程为%（2）⽤拉格朗⽇法插值%定义Lagrange程序function f=Language(x,y,x0)syms t;if(length(x)==length(y))n=length(x);elsedisp('xoíyµêy2??àµè￡?');return;endf=0.0;for(i=1:n)l=y(i);for(j=1:i-1)l=l*(t-x(j))/(x(i)-x(j));end ;for (j=i+1:n)l=l*(t-x(j))/(x(i)-x(j)); end ; f=f+l; simplify(f); if (i==n)if (nargin==3) f=subs(f,'t',x0); elsef=collect(f); f=vpa(f,6); end end end>> Language(x,y) ans =52462.6*t+189995.*t^3-189851.*t^4+136778.*t^5-11.3161*t^12-.277283e-6*t^18+1.18284*t^13-73866.6*t^6+.111076e-4*t^17-.976904e-1*t^14+.427949e-8*t^19-.307453e-10*t^20+30677.6*t^7+2564.20*t^9-9968.98*t^8+.628590e-2*t^15-525.813*t^10-9652.78-.308159e-3*t^16+86.2514*t^11-128683.*t^2⼆、已知Wilson 矩阵1078775658610975910A=，且向量32233331b ??=，则⽅程组Ax b =有准确解[]1111Tx =。
数值分析实验2014
数值分析实验(2014,9,16~10,28)信计1201班，人数34人数学系机房数值分析计算实习报告册专业__________________学号_______________姓名_______________2014~2015年第一学期实验一数值计算的工具Matlab1. 解释下MATLABS序的输出结果程序：t=0.1n=1:10e=n/10-n*te 的结果：0 0 -5.5511e-017 0 0-1.1102e-016 -1.1102e-016 0 0 02. 下面MATLABS序的的功能是什么？程序：x=1;while 1+x>1,x=x/2,pause(0.02),e nd用迭代法求出x=x/2,的最小值x=1;while x+x>x,x=2*x,pause(0.02),e nd用迭代法求出x=2*x,的值，使得2x>Xx=1;while x+x>x,x=x/2,pause(0.02),e nd用迭代法求出x=x/2,的最小值，使得2x>X3. 考虑下面二次代数方程的求解问题2ax bx c = 0公式x=电上4ac是熟知的，与之等价地有_____________________________ ,对于2a-b ■ b -4aca =1,b =100000000,c =1，应当如何选择算法。
b ~4ac计算，因为b与b2— 4ac相近，两个相加减不宜应该用2a u做分母3 5 74. 函数sin(x)有幂级数展开sin x = x - x - - ■■3! 5! 7!利用幕级数计算sinx的MATLAB程序为fun cti on s=powers in(x)s=0;t=x;n=1;while s+t~=s;s=s+t ；t=-x A2/ ((n+1)*(n+2) ) *t ;n=n+2 ；endt仁cputime;pause(10);t2=cputime;t0=t2-t1(a) 解释上述程序的终止准则。
数值分析计算实习题
插值法1.下列数据点的插值x 0 1 4 9 16 25 36 49 64y 0 1 2 3 4 5 6 7 8可以得到平方根函数的近似，在区间[0,64] 上作图.(1)用这9 个点作8 次多项式插值Ls(x).(2)用三次样条( 第一边界条件)程序求S(x).从得到结果看在[0,64] 上，哪个插值更精确；在区间[0,1] 上，两种插值哪个更精确解：(1) 拉格朗日插值多项式，求解程序如下syms x l;x1=[0 1 4 9 16 25 36 49 64]; y1=[0 1 2 3 4 5 6 7 8]; n=length(x1);Ls=sym(0);for i=1:nl=sym(y1(i));for k=1:i-1l=l*(x-x1(k))/(x1(i)-x1(k));endfor k=i+1:n l=l*(x-x1(k))/(x1(i)-x1(k));endLs=Ls+l; endLs=simplify(Ls) % 为所求插值多项式Ls(x).输出结果为Ls = -/*xA2+95549/72072*x-1/00*xA8-2168879/0*xA4+19/0*xA7+657859/*xA3+33983/ 0*xA5-13003/00*xA6(2) 三次样条插值，程序如下x1=[0 1 4 9 16 25 36 49 64];y1=[0 1 2 3 4 5 6 7 8];x2=[0:1:64];y3=s plin e(x1,y1,x2);p=po Iyfit(x2,y3,3); % 得到三次样条拟合函数S=p(1)+p(2)*x+ p(3)*x^2+p(4)*xA3 % 得到S(x) 输出结果为：S =/6464-2399/88*x+/1984*xA2+2656867/624*xA3⑶ 在区间［0,64］上，分别对这两种插值和标准函数作图,Plot(x2,sqrt(x2),'b',x2,y2,'r',x2,y3,'y')蓝色曲线为y="X函数曲线，红色曲线为拉格朗日插值函数曲线，黄色曲线为三次样条插值曲线可以看到蓝色曲线与黄色曲线几乎重合，因此在区间[0,64] 上三次样条插值更精确。
数学与应用数学实习生的实习总结
数学与应用数学实习生的实习总结在过去的几个月里，作为一名数学与应用数学专业的实习生，我有幸在一家著名科研机构进行了为期三个月的实习。
在这段时间内，我不仅学到了很多数学理论知识，还获得了宝贵的实践经验。
本文将围绕我在实习期间所参与的项目和我所取得的成果进行总结，并分享一些在实习过程中的收获和感悟。
首先，在实习的第一个月里，我主要参与了一个关于图论的研究项目。
这个项目旨在研究图的最大团问题，并找到一种高效的算法来解决这个问题。
在导师的指导下，我深入学习了图论的基本概念和算法，并利用Python编程语言实现了一个基于回溯法的求解最大团问题的程序。
我通过运行实例，验证了该算法的正确性和高效性。
这个项目不仅让我对图论有了更深入的了解，还提升了我的编程能力和问题解决能力。
其次，在实习的第二个月里，我加入了一个数值计算的研究小组。
我们小组的研究方向是使用有限元方法求解偏微分方程的数值解。
在这个项目中，我学习了有限元方法的基本原理和常用算法，并使用MATLAB编程语言实现了一个二维热传导方程的有限元求解程序。
通过对程序的不断调试和优化，我成功求解了一些实际问题，并将结果与解析解进行了比较，验证了有限元方法在求解偏微分方程数值解中的可靠性和有效性。
通过参与这个项目，我对数值计算领域有了更深入的认识，并锻炼了我的编程和数学建模能力。
另外，在实习的最后一个月里，我参与了一个关于金融风险控制的研究项目。
这个项目的目标是建立一个风险评估模型，以帮助金融机构更好地管理风险。
在这个项目中，我学习了一些金融风险管理的基本理论和方法，包括价值-at-风险（VaR）和条件价值-at-风险（CVaR）等。
我利用R编程语言实现了这些风险评估模型，并使用实际数据进行了模拟实验和风险评估。
通过这个项目，我对金融风险控制有了更深入的理解，同时也提升了我的数据处理和分析能力。
在整个实习期间，我不仅仅学到了专业知识和技能，还获得了一些宝贵的经验和感悟。
数值分析第五版计算实习题
弟二草插值法3.卜列数据点的插值可以得到平方根函数的近似，在区间064］上作图。
(1〉用这9个点做8次多项式插值Q x)。
(2)用三次样条(第一边界条件)程岸求S(X)。
从得到结果石在［0.64］ 1：・哪个插值更粘确：在区间［0,1］ I：•两种插值哪个更精确?(1) 8次多项式插值：(1)8次多项式插值:首先建立新的M-file:输入如卜代码(此为拉格朗口插值的功能函数)并保存function f=Language(x,y,x0)%求Li知数据点的拉格朗Fl插值多项式%己知数据点的x坐标向量：x%已知数据点的y坐标向量:y%插值的x坐标：x0%求得的拉格朗H插值多项式或在X0处的插值：fsyms t;ifi(lcngth(x)=length(y))n=length(x);elsedisp(*x和y的维数不相等!)；return;end %检错tbr(i=l:n)i=y(i);fbr(j=1:i-l)l=l*(t-x(j))/(x(i)-x(j))；end;for(j=i-M:n)end;for(j=i+l:n) l=l*(t-x(j))/(x(i)-x(j))； end;simplify(f);if(i==n) if|nargin=3)f=subs(C't\xO);else f=collcct(f);f=vpa(f,6);endendend再建立新的M-file：输入：clear;x=[0 1 49 16 25 36 49 64];y=[0：l：8]；%计算拉格朗口基丞数%计算拉格朗ri插值函数%化简%计算插值点的曲数值%将插值多项式展开%将插值多项式的系数化成6位精度的小数f=Uinguage(x,y) 运行得到f=1.32574*1-381410*t A2+.604294e-1 *t A3+.222972e-3 *t A5-.542921 e-5*t A6+.671268e・7T7・.328063e・9T8・.498071 e-2*t A4 这就是8次多项式插值L s(x)= 1.32574怜.381410*t A2+.604294e-1 *t A3+.222972e-3 *t A5-.542921 e-5*t A6+.671268e-7*t A7-.328063e-9*t A8-.498071 e-2*t A4. (2)三次样条插值：建立新的M-filc：输入：clear;x=[0 I 49 1625 36 4964];尸[0:8];t=[0:0.1:64];Y=t.A(0.5);O=Language(x,y)f= 1,32574*t-.381410*t.A2+.604294e-1 *t.A3+.222972e-3*t.A5-.542921 e・5*(. W+.671268e-7*t.A7-.328063e-9*t.A8-.498071 e-2 *t.A4;S=interp l(x,y,t.'spline,);plol(x,y,o；(・YY.lf.'b'」S'g:');grid;运行程序得到如下图：从结果屮很明显可以看出在[0.64].上.三次样条插值更精确，儿乎与原函数帀合。
数值分析计算实习题答案
数值分析计算实习题答案数值分析计算实习题答案数值分析是一门研究如何利用计算机对数学问题进行近似求解的学科。
在数值分析的学习过程中，实习题是一种重要的学习方式，通过实践来巩固理论知识，并培养解决实际问题的能力。
本文将为大家提供一些数值分析计算实习题的答案，希望能够帮助大家更好地理解和掌握数值分析的相关知识。
一、插值与拟合1. 已知一组数据点，要求通过这些数据点构造一个一次插值多项式，并求出在某一特定点的函数值。
答案：首先，我们可以根据给定的数据点构造一个一次插值多项式。
假设给定的数据点为(x0, y0), (x1, y1)，我们可以构造一个一次多项式p(x) = a0 + a1x，其中a0和a1为待定系数。
根据插值条件，我们有p(x0) = y0，p(x1) = y1。
将这两个条件代入多项式中，可以得到一个方程组，通过求解这个方程组，我们就可以确定a0和a1的值。
最后，将求得的多项式代入到某一特定点，就可以得到该点的函数值。
2. 已知一组数据点，要求通过这些数据点进行最小二乘拟合，并求出拟合曲线的表达式。
答案：最小二乘拟合是一种通过最小化误差平方和来找到最佳拟合曲线的方法。
假设给定的数据点为(x0, y0), (x1, y1)，我们可以构造一个拟合曲线的表达式y =a0 + a1x + a2x^2 + ... + anx^n，其中a0, a1, ..., an为待定系数。
根据最小二乘拟合原理，我们需要最小化误差平方和E = Σ(yi - f(xi))^2，其中yi为实际数据点的y值，f(xi)为拟合曲线在xi处的函数值。
通过求解这个最小化问题，我们就可以确定拟合曲线的表达式。
二、数值积分1. 已知一个函数的表达式，要求通过数值积分的方法计算函数在某一区间上的定积分值。
答案：数值积分是一种通过将定积分转化为数值求和来近似计算的方法。
假设给定的函数表达式为f(x)，我们可以将定积分∫[a, b]f(x)dx近似为Σwi * f(xi)，其中wi为权重系数，xi为待定节点。
数值分析实习报告总结
一、实习背景数值分析是数学的一个重要分支，它研究如何用数值方法求解数学问题。
随着计算机技术的飞速发展，数值分析在各个领域得到了广泛的应用。
为了提高自己的实践能力，我选择了数值分析作为实习课题，希望通过这次实习，能够掌握数值分析的基本方法，并将其应用于实际问题中。
二、实习过程1. 实习初期在实习初期，我首先了解了数值分析的基本概念、理论和方法。
通过阅读相关教材和文献，我对数值分析有了初步的认识。
接着，我学习了数值分析的基本方法，如泰勒展开、牛顿法、高斯消元法等。
2. 实习中期在实习中期，我选择了几个实际问题进行数值计算。
首先，我使用泰勒展开法求解一个简单的微分方程。
通过编写程序，我得到了微分方程的近似解。
然后，我运用牛顿法求解一个非线性方程组。
在实际计算过程中，我遇到了一些问题，如收敛性、迭代次数过多等。
通过查阅资料和请教导师，我找到了解决方法，成功求解了方程组。
3. 实习后期在实习后期，我进一步学习了数值分析的高级方法，如复化梯形公式、复化Simpson公式、自适应梯形法等。
这些方法在解决实际问题中具有更高的精度和效率。
我选择了一个具体的工程问题，运用复化梯形公式求解定积分。
在计算过程中，我遇到了区间细分、精度控制等问题。
通过不断尝试和调整，我得到了较为精确的积分值。
三、实习收获与体会1. 理论与实践相结合通过这次实习，我深刻体会到理论与实践相结合的重要性。
在实习过程中，我不仅学习了数值分析的理论知识，还将其应用于实际问题中。
这使我更加深刻地理解了数值分析的基本方法，提高了自己的实践能力。
2. 严谨的学术态度在实习过程中，我养成了严谨的学术态度。
在编写程序、进行数值计算时，我注重细节，力求精确。
这使我更加注重学术规范，提高了自己的学术素养。
3. 团队合作精神实习过程中，我与其他同学进行了交流与合作。
在解决实际问题时，我们互相学习、互相帮助，共同完成了实习任务。
这使我更加懂得团队合作的重要性，提高了自己的团队协作能力。
北航数值分析计算实习题目一 幂法 反幂法 求矩阵特征值
《数值分析》计算实习题目第一题：1. 算法设计方案（1）1λ，501λ和s λ的值。
1)首先通过幂法求出按模最大的特征值λt1，然后根据λt1进行原点平移求出另一特征值λt2，比较两值大小，数值小的为所求最小特征值λ1，数值大的为是所求最大特征值λ501。
2)使用反幂法求λs ，其中需要解线性方程组。
因为A 为带状线性方程组，此处采用LU 分解法解带状方程组。
（2）与140k λλμλ-5011=+k 最接近的特征值λik 。
通过带有原点平移的反幂法求出与数k μ最接近的特征值 λik 。
（3）2cond(A)和det A 。
1）1=nλλ2cond(A)，其中1λ和n λ分别是按模最大和最小特征值。
2）利用步骤（1）中分解矩阵A 得出的LU 矩阵，L 为单位下三角阵,U 为上三角阵，其中U 矩阵的主对角线元素之积即为det A 。
由于A 的元素零元素较多，为节省储存量，将A 的元素存为6×501的数组中，程序中采用get_an_element()函数来从小数组中取出A 中的元素。
2.全部源程序#include <stdio.h>#include <math.h>void init_a();//初始化Adouble get_an_element(int,int);//取A 中的元素函数double powermethod(double);//原点平移的幂法double inversepowermethod(double);//原点平移的反幂法int presolve(double);//三角LU 分解int solve(double [],double []);//解方程组int max(int,int);int min(int,int);double (*u)[502]=new double[502][502];//上三角U 数组double (*l)[502]=new double[502][502];//单位下三角L 数组double a[6][502];//矩阵Aint main(){int i,k;double lambdat1,lambdat2,lambda1,lambda501,lambdas,mu[40],det;double lambda[40];init_a();//初始化Alambdat1=powermethod(0);lambdat2=powermethod(lambdat1);lambda1=lambdat1<lambdat2?lambdat1:lambdat2;lambda501=lambdat1>lambdat2?lambdat1:lambdat2;presolve(0);lambdas=inversepowermethod(0);det=1;for(i=1;i<=501;i++)det=det*u[i][i];for (k=1;k<=39;k++){mu[k]=lambda1+k*(lambda501-lambda1)/40;presolve(mu[k]);lambda[k]=inversepowermethod(mu[k]);}printf("------------所有特征值如下------------\n");printf("λ=%1.11e λ=%1.11e\n",lambda1,lambda501);printf("λs=%1.11e\n",lambdas);printf("cond(A)=%1.11e\n",fabs(lambdat1/lambdas));printf("detA=%1.11e \n",det);for (k=1;k<=39;k++){printf("λi%d=%1.11e ",k,lambda[k]);if(k % 3==0) printf("\n");} delete []u;delete []l;//释放堆内存return 0;}void init_a()//初始化A{int i;for (i=3;i<=501;i++) a[1][i]=a[5][502-i]=-0.064;for (i=2;i<=501;i++) a[2][i]=a[4][502-i]=0.16;for (i=1;i<=501;i++) a[3][i]=(1.64-0.024*i)*sin(0.2*i)-0.64*exp(0.1/i); }double get_an_element(int i,int j)//从A中节省存储量的提取元素方法{if (fabs(i-j)<=2) return a[i-j+3][j];else return 0;}double powermethod(double offset)//幂法{int i,x1;double u[502],y[502];double beta=0,prebeta=-1000,yita=0;for (i=1;i<=501;i++)u[i]=1,y[i]=0;//设置初始向量u[]for (int k=1;k<=10000;k++){yita=0;for (i=1;i<=501;i++) yita=sqrt(yita*yita+u[i]*u[i]);for (i=1;i<=501;i++) y[i]=u[i]/yita;for (x1=1;x1<=501;x1++){u[x1]=0;for (int x2=1;x2<=501;x2++)u[x1]=u[x1]+((x1==x2)?(get_an_element(x1,x2)-offset):get_an_element(x1,x2))*y[x2] ;}prebeta=beta;beta=0;for (i=1;i<=501;i++) beta=beta+ y[i]*u[i];if (fabs((prebeta-beta)/beta)<=1e-12) {printf("offset=%f lambda=%f err=%e k=%d\n",offset,(beta+offset),fabs((prebeta-beta)/beta),k);break;};//输出中间过程,包括偏移量,误差,迭代次数}return (beta+offset);}double inversepowermethod(double offset)//反幂法{int i;double u[502],y[502];double beta=0,prebeta=0,yita=0;for (i=1;i<=501;i++)u[i]=1,y[i]=0; //设置初始向量u[]for (int k=1;k<=10000;k++){yita=0;for (i=1;i<=501;i++) yita=sqrt(yita*yita+u[i]*u[i]);for (i=1;i<=501;i++) y[i]=u[i]/yita;solve(u,y);prebeta=beta;beta=0;for (i=1;i<=501;i++) beta=beta+ y[i]*u[i];beta=1/beta;if (fabs((prebeta-beta)/beta)<=1e-12) {printf("offset=%f lambda=%f err=%e k=%d\n",offset,(beta+offset),fabs((prebeta-beta)/beta),k);break;};//输出中间过程,包括偏移量,误差,迭代次数}return (beta+offset);int presolve(double offset)//三角LU分解{int i,k,j,t;double sum;for (k=1;k<=501;k++)for (j=1;j<=501;j++){u[k][j]=l[k][j]=0;if (k==j) l[k][j]=1;} //初始化LU矩阵for (k=1;k<=501;k++){for (j=k;j<=min(k+2,501);j++){sum=0;for (t=max(1,max(k-2,j-2)) ; t<=(k-1) ; t++)sum=sum+l[k][t]*u[t][j];u[k][j]=((k==j)?(get_an_element(k,j)-offset):get_an_element(k,j))-sum;}if (k==501) continue;for (i=k+1;i<=min(k+2,501);i++){sum=0;for (t=max(1,max(i-2,k-2));t<=(k-1);t++)sum=sum+l[i][t]*u[t][k];l[i][k]=(((i==k)?(get_an_element(i,k)-offset):get_an_element(i,k))-sum)/u[k][k];}}return 0;}int solve(double x[],double b[])//解方程组{int i,t;double y[502];double sum;y[1]=b[1];for (i=2;i<=501;i++){sum=0;for (t=max(1,i-2);t<=i-1;t++)sum=sum+l[i][t]*y[t];y[i]=b[i]-sum;}x[501]=y[501]/u[501][501];for (i=500;i>=1;i--){sum=0;for (t=i+1;t<=min(i+2,501);t++)sum=sum+u[i][t]*x[t];x[i]=(y[i]-sum)/u[i][i];}return 0;}int max(int x,int y){return (x>y?x:y);}int min(int x,int y){return (x<y?x:y);}3.计算结果结果如下图所示:部分中间结果：给出了偏移量(offset)，误差(err)，迭代次数(k)4.讨论迭代初始向量的选取对计算结果的影响,并说明原因使用u[i]=1（i=1,2,...,501）作为初始向量进行迭代，可得出以上结果。
数值分析计算实习第一题
直接用定义: ������������(������������)2 = ‖������������‖2‖������������−1‖2
求 A 的条件数很繁琐,需要先进行化简：
首先：
由于 A 是对称矩阵，
‖������������‖2 = �������������max(������������������������������������)
说明 ：
1. 在所用的算法中，凡是要给出精度水平的ε，都取 ������������=10−12。
2. 选择算法的时候应使矩阵 A 的所有零元素都不存储。
3. 打印以下内容：
(1)算法设计方案和思路。
(2)全部源程序。
(3)特征值������������1,������������501,������������������������,������������������������������������(������������=1,2,⋯,39)以及������������������������������������������������(������������)2, det������������的值(采用 e 型输出实型数，并 至少显示 12 位有效数字)。
λi[16] -2.533970311130E+00 λi[38] 8.648666065193E+00
λi[17] -2.003230769563E+00 λi[39] 9.254200344575E+00
λi[18] -1.503557611227E+00 cond(A)2 1.925204273903E+03
λi[19] -9.935586060080E-01 det(A) 2.772786141752E+118
数值分析报告(第五版)计算实习的题目第三章
数值分析计算实习题第三章第二次作业：题一：x=-1:0.2:1;y=1./(1+25.*x.^2);f1=polyfit(x,y,3)f=poly2sym(f1)y1=polyval(f1,x)x2=linspace(-1,1,10)y2=interp1(x,y,x2)plot(x,y,'r*-',x,y1,'b-')hold onplot(x2,y2,'k')legend('数据点','3次拟合曲线','3次多项式插值')xlabel('X'),ylabel('Y')输出：f1 =0.0000 -0.5752 0.0000 0.4841f =(4591875547102675*x^3)/81129638414606681695789005144064 - (3305*x^2)/5746 + (1469057404776431*x)/20282409603651670423947251286016 + 4360609662300613/9007199254740992y1 =-0.0911 0.1160 0.2771 0.3921 0.4611 0.4841 0.4611 0.3921 0.2771 0.1160 -0.0911x2 =-1.0000 -0.7778 -0.5556 -0.3333 -0.1111 0.1111 0.3333 0.5556 0.7778 1.0000y2 =0.0385 0.0634 0.1222 0.3000 0.7222 0.7222 0.3000 0.1222 0.0634 0.0385题二：X=[0.0 0.1 0.2 0.3 0.5 0.8 1.0];Y=[1.0 0.41 0.50 0.61 0.91 2.02 2.46];p1=polyfit(X,Y,3)p2=polyfit(X,Y,4)Y1=polyval(p1,X)Y2=polyval(p2,X)plot(X,Y,'r*',X,Y1,'b-.',X,Y2,'g--')p3=polyfit(X,Y,2)Y3=polyval(p3,X)f1=poly2sym(p1)f2=poly2sym(p2)f3=poly2sym(p3)plot(X,Y,'r*',X,Y1,'b-.',X,Y2,'g--',X,Y3,'m--')legend('数据点','3次多项式拟合','4次多项式拟合','2次多项式拟合') xlabel('X轴'),ylabel('Y轴')输出：p1 =-6.6221 12.8147 -4.6591 0.9266p2 =2.8853 -12.3348 16.2747 -5.2987 0.9427Y1 =0.9266 0.5822 0.4544 0.5034 0.9730 2.0103 2.4602Y2 =0.9427 0.5635 0.4399 0.5082 1.0005 1.9860 2.4692p3 =3.1316 -1.2400 0.7356Y3 =0.7356 0.6429 0.6128 0.6454 0.8984 1.7477 2.6271f1 =- (7455778416425075*x^3)/1125899906842624 + (1803512222945435*x^2)/140737488355328 - (40981580032809*x)/8796093022208 + 8345953784399011/9007199254740992f2 =(1624271450198125*x^4)/562949953421312 - (3471944732519173*x^3)/281474976710656 + (4580931990070659*x^2)/281474976710656 - (1491459232922115*x)/281474976710656 + 1061409433081293/1125899906842624f3 =(18733*x^2)/5982 - (74179*x)/59820 + 73337/99700题三：建立三角插值函数的m文件function [A,B,Y1,Rm]=sanjiaobijin(X,Y,X1,m)%A B分别是m阶三角多项式Tm （x）的系数aj，bj（j=1,2,...,m）的系数矩阵，Y1是Tm（x）在X1处的值，X Y 数据点 ,Rm为均方误差n=length(X)-1;max1=fix((n-1)/2);if m>max1m=max1;endA=zeros(1,m+1);B=zeros(1,m+1);Ym=(Y(1)+Y(n+1))/2;Y(1)=Ym;Y(n+1)=Ym;A(1)=2*sum(Y)/n;for i=1:mB(i+1)=sin(i*X)*Y';A(i+1)=cos(i*X)*Y';endA=2*A/n;B=2*B/n;A(1)=A(1)/2;Y1=A(1);for k=1:mY1=Y1+A(k+1)*cos(k*X1)+B(k+1)*sin(k*X1);Tm=A(1)+A(k+1).*cos(k*X)+B(k+1).*sin(k*X);k=k+1;endY,Tm,Rm=(sum(Y-Tm).^2)/n输出：>> X=-pi:2*pi/33:pi;>> Y=X.^2.*cos(X);[A,B,Y1,Rm]=sanjiaobijin(X,Y,X1,16)输出：A =1 至 12 列-0.1397 4.4002 -2.8326 1.2355 -0.9128 0.7914 -0.7319 0.6982 -0.6773 0.6635 -0.6541 0.647413 至 17 列-0.6426 0.6393 -0.6370 0.6355 -0.6348B =1.0e-15 *1 至 12 列0 -0.0194 -0.0150 -0.0044 -0.0300 0.0105 0.0627 -0.0821 -0.0599 -0.0133 -0.0211 0.029713 至 17 列0.0178 0.0962 -0.1049 0.0328 -0.0122即可得16插值多项式的值X1=-pi:0.001:pi;[A,B,Y1,Rm]=sanjiaobijin(X,Y,X1,16)plot(X,Y,'r*',X1,Y1,'b-.')legend('数据点','16次三角插值多项式')xlabel('X轴'),ylabel('Y轴')。
数值计算方法上机实习题答案.doc
1．设I n 1 x ndx ，0 5 x（ 1）由递推公式 I n 5I n 11，从 I 0的几个近似值出发，计算I 20；n解：易得： I 0 ln6-ln5=0.1823, 程序为：I=0.182;for n=1:20I=(-5)*I+1/n;endI输出结果为： I 20= -3.0666e+010（ 2）粗糙估计 I 20，用 I n 1 1I n 1 1 ，计算 I 0；5 5n0.0079 1 x 20 1 x 200.0095因为dx I 20dx 6 5所以取 I 20 1(0.0079 0.0095) 0.0087 2程序为： I=0.0087;for n=1:20I=(-1/5)*I+1/(5*n);endII 0= 0.0083（ 3）分析结果的可靠性及产生此现象的原因(重点分析原因 )。
首先分析两种递推式的误差；设第一递推式中开始时的误差为E0 I 0 I 0，递推过程的舍入误差不计。
并记 E n I n I n，则有 E n 5E n 1 ( 5) n E0。
因为 E20( 5) 20 E0 I 20，所此递推式不可靠。
而在第二种递推式中E0 1E1 (1)n E n，误差在缩小，5 5所以此递推式是可靠的。
出现以上运行结果的主要原因是在构造递推式过程中，考虑误差是否得到控制，即算法是否数值稳定。
2．求方程e x10x 2 0 的近似根，要求x k 1x k 5 10 4，并比较计算量。
（1）在 [0， 1]上用二分法；程序： a=0;b=1.0;while abs(b-a)>5*1e-4c=(b+a)/2;if exp(c)+10*c-2>0b=c;else a=c;endendc结果： c =0.0903（ 2）取初值x0 0，并用迭代 x k 1 2 e x ；10程序： x=0;a=1;while abs(x-a)>5*1e-4a=x;x=(2-exp(x))/10;endx结果： x =0.0905（3）加速迭代的结果；程序： x=0;a=0;b=1;while abs(b-a)>5*1e-4a=x;y=exp(x)+10*x-2;z=exp(y)+10*y-2;x=x-(y-x)^2/(z-2*y+x);b=x;endx结果： x =0.0995（ 4）取初值x00 ，并用牛顿迭代法；程序： x=0;a=0;b=1;while abs(b-a)>5*1e-4a=x;x=x-(exp(x)+10*x-2)/(exp(x)+10); b=x;end x 结果： x =0.0905（ 5） 分析绝对误差。
科学计算软件开发：实习中的数值模拟与算法实现
科学计算软件开发：实习中的数值模拟与算法实现科学计算软件开发是现代科技领域不可或缺的重要组成部分，它广泛应用于物理学、化学、工程学等领域的数值模拟和算法实现。
在实习中，我有幸参与了一项科学计算软件的开发工作，亲身体验了数值模拟与算法实现的具体过程。
本文将重点介绍我在实习中所参与的数值模拟项目，并分享一些经验与感悟。
在实习中，我们的项目是基于Python语言开发的数值模拟软件。
这个软件是用于模拟光传输过程的，基于有限差分方法和蒙特卡洛方法构建了一个精细的光传输模型，并利用算法对光在不同介质中的传递、反射、折射等过程进行模拟。
我们的目标是为科学研究和工程应用提供一个可靠、高效的数值模拟工具。
首先，我们需要对模拟问题进行建模。
在我们的项目中，光的传输过程可以用波动方程进行描述。
我们将空间离散化为小立方体，时间离散化为小时间步长。
通过有限差分方法，我们可以将波动方程离散化为一组线性方程，每个方程表示一个立方体上的光强度。
这样，我们就得到了一个庞大的线性方程组，通过求解这个方程组，就能够得到光的传输过程的数值解。
而蒙特卡洛方法则用于模拟光的随机传输过程，通过随机采样和统计平均的方法得到光的传输特性。
接下来，我们需要实现这些算法。
对于有限差分方法，我们使用Python的科学计算库NumPy来进行矩阵运算，通过构造系数矩阵和右端向量，利用线性代数库解线性方程组。
对于蒙特卡洛方法，我们需要生成服从特定概率分布的随机数，并对大量的随机样本进行模拟。
为了提高计算效率，我们使用了Python的并行计算库multiprocessing来并行计算多个样本。
在实习中，我主要负责算法的实现和结果的分析。
在实现算法时，我遇到了很多挑战和困惑。
首先，由于我是计算机专业的学生，对于光传输过程和相关物理知识并不是很熟悉，所以需要花费很多时间来学习相关的知识。
其次，由于算法的复杂性和规模较大，我需要对代码进行优化和调试，以提高计算效率和准确性。
数值分析计算实习题
1.下列数据点的插值x 0 1 4 9 16 25 36 49 64y 0 1 2 3 4 5 6 7 8可以得到平方根函数的近似，在区间[0,64]上作图.(1)用这9个点作8次多项式插值Ls(x).(2)用三次样条(第一边界条件)程序求S(x).从得到结果看在[0,64]上，哪个插值更精确；在区间[0,1]上，两种插值哪个更精确？解：(1)拉格朗日插值多项式，求解程序如下syms x l;x1=[0 1 4 9 16 25 36 49 64];y1=[0 1 2 3 4 5 6 7 8];n=length(x1);Ls=sym(0);for i=1:nl=sym(y1(i));for k=1:i-1l=l*(x-x1(k))/(x1(i)-x1(k));endfor k=i+1:nl=l*(x-x1(k))/(x1(i)-x1(k));endLs=Ls+l;endLs=simplify(Ls) %为所求插值多项式Ls(x).输出结果为Ls =-/*x^2+95549/72072*x-1/00*x^8-2168879/0*x^4+19/0*x^7+657859/*x^3+33983/ 0*x^5-13003/00*x^6(2)三次样条插值，程序如下x1=[0 1 4 9 16 25 36 49 64];y1=[0 1 2 3 4 5 6 7 8];x2=[0:1:64];y3=spline(x1,y1,x2);p=polyfit(x2,y3,3); %得到三次样条拟合函数S=p(1)+p(2)*x+p(3)*x^2+p(4)*x^3 %得到S(x)输出结果为：S =/6464-2399/88*x+/1984*x^2+2656867/624*x^3(3)在区间[0,64]上，分别对这两种插值和标准函数作图，plot(x2,sqrt(x2),'b',x2,y2,'r',x2,y3,'y')蓝色曲线为y=函数曲线,红色曲线为拉格朗日插值函数曲线，黄色曲线为三次样条插值曲线可以看到蓝色曲线与黄色曲线几乎重合，因此在区间[0,64]上三次样条插值更精确。
数值分析第五版计算实习题
数值分析计算实习题第二章2-1程序：clear;clc;x1=[0.2 0.4 0.6 0.8 1.0];y1=[0.98 0.92 0.81 0.64 0.38];n=length(y1);c=y1(:);for j=2:n %求差商for i=n:-1:jc(i)=(c(i)-c(i-1))/(x1(i)-x1(i-j+1));endendsyms x df d;df(1)=1;d(1)=y1(1);for i=2:n %求牛顿差值多项式df(i)=df(i-1)*(x-x1(i-1));d(i)=c(i)*df(i);enddisp('4次牛顿插值多项式');P4=vpa(collect((sum(d))),5) %P4即为4次牛顿插值多项式,并保留小数点后5位数pp=csape(x1,y1, 'variational');%调用三次样条函数q=pp.coefs;disp('三次样条函数');for i=1:4S=q(i,:)*[(x-x1(i))^3;(x-x1(i))^2;(x-x1(i));1];S=vpa(collect(S),5)endx2=0.2:0.08:1.08;dot=[1 2 11 12];figureezplot(P4,[0.2,1.08]);hold ony2=fnval(pp,x2);x=x2(dot);y3=eval(P4);y4=fnval(pp,x2(dot));plot(x2,y2,'r',x2(dot),y3,'b*',x2(dot),y4,'co');title('4次牛顿插值及三次样条');结果如下：4次牛顿插值多项式P4 = - 0.52083*x^4 + 0.83333*x^3 - 1.1042*x^2 + 0.19167*x + 0.98三次样条函数x∈[0.2,0.4]时， S = - 1.3393*x^3 + 0.80357*x^2 - 0.40714*x + 1.04 x∈[0.4,0.6]时，S = 0.44643*x^3 - 1.3393*x^2 + 0.45*x + 0.92571x∈[0.6,0.8]时，S = - 1.6964*x^3 + 2.5179*x^2 - 1.8643*x + 1.3886 x∈[0.8,1.0]时，S =2.5893*x^3 - 7.7679*x^2 + 6.3643*x - 0.80571输出图如下2-3（1）程序：clear;clc;x1=[0 1 4 9 16 25 36 49 64];y1=[0 1 2 3 4 5 6 7 8];%插值点n=length(y1);a=ones(n,2);a(:,2)=-x1';c=1;for i=1:nc=conv(c,a(i,:));endq=zeros(n,n);r=zeros(n,n+1);for i=1:n[q(i,:),r(i,:)]=deconv(c,a(i,:));%wn+1/(x-xk) end 三次样条插值曲线4次牛顿插值曲线Dw=zeros(1,n);for i=1:nDw(i)=y1(i)/polyval(q(i,:),x1(i));%系数endp=Dw*q;syms x L8;for i=1:nL8(i)=p(n-i+1)*x^(i-1);enddisp('8次拉格朗日插值');L8=vpa(collect((sum(L8))),5)xi=0:64;yi=polyval(p,xi);figureplot(xi,yi,x1,y1,'r*');hold ontitle('8次拉格朗日插值');结果如下：8次拉格朗日插值L8 =- 3.2806e-10*x^8 + 6.7127e-8*x^7 - 5.4292e-6*x^6 + 0.00022297*x^5 - 0.0049807*x^4 + 0.060429*x^3 - 0.38141*x^2 + 1.3257*x输出图如下：第五章4-1（3）程序：clc;clear;y= @(x) sqrt(x).*log(x);a=0;b=1;tol=1e-4;p=quad(y,a,b,tol);fprintf('采用自适应辛普森积分结果为： %d \n', p);结果如下：采用自适应辛普森积分结果为： -4.439756e-01第九章9-1(a)程序：clc;clear;a=1;b=2;%定义域h=0.05;%步长n=(b-a)/h;y0=1;%初值f= @(x,y) 1/x^2-y/x;%微分函数Xn=linspace(a,b,n+1);%将定义域分为n等份Yn=zeros(1,n);%结果矩阵Yn(1)=y0;%赋初值%以下根据改进欧拉公式求解for i=1:nxn=Xn(i);xnn=Xn(i+1);yn=Yn(i);yp=yn+h*f(xn,yn);yc=yn+h*f(xnn,yp);yn=(yp+yc)/2;Yn(i+1)=yn;endXn=Yn;%以下根据经典四阶R-K法公式求解for i=1:nxn=Xn(i);yn=Yn(i);k1=f(xn,yn);k2=f(xn+h/2,yn+h/2*k1);k3=f(xn+h/2,yn+h/2*k2);k4=f(xn+h,yn+h*k3);yn=yn+h/6*(k1+2*k2+2*k3+k4);Yn(i+1)=yn;enddisp(' 改进欧拉法四阶经典R-K法'); disp([Xn' Yn'])结果如下：改进欧拉法四阶经典R-K法1 10.99887 0.998850.99577 0.99780.99114 0.996940.98532 0.996340.97857 0.996030.97111 0.996060.96311 0.996450.9547 0.997230.94598 0.998410.93705 10.92798 1.0020.91883 1.00440.90964 1.00730.90045 1.01060.89129 1.01430.88218 1.01840.87315 1.02290.86421 1.02780.85538 1.03310.84665 1.0388（b）程序：clc;clear;a=0;b=1;%定义域H=[0.1 0.025 0.01];%步长y0=1/3;%初值f= @(x,y) -50*y+50*x^2+2*x;%微分函数xi=linspace(a,b,11);Y=1/3*exp(-50*xi)+xi.^2;%准确解Ym=zeros(1,11);for j=1:3h=H(j);n=(b-a)/h;Xn=linspace(a,b,n+1);%将定义域分为n等份Yn=zeros(1,n);%结果矩阵Yn(1)=y0;%赋初值for i=1:nxn=Xn(i);yn=Yn(i);k1=f(xn,yn);k2=f(xn+h/2,yn+h/2*k1);k3=f(xn+h/2,yn+h/2*k2);k4=f(xn+h,yn+h*k3);yn=yn+h/6*(k1+2*k2+2*k3+k4);Yn(i+1)=yn;endfor k=1:11m=0.1/h;Ym(k)=Yn(1+(k-1)*m);enddelta=Ym-Y;fprintf('步长为： %d \n', h);disp(' 四阶经典R-K法准确解误差'); disp([Ym' Y' delta'])end结果如下：步长为： 1.000000e-01四阶经典R-K法准确解误差0.33333 0.33333 04.6055 0.012246 4.593263.062 0.040015 63.022864.05 0.09 863.9611844 0.16 118431.6235e+05 0.25 1.6235e+052.2256e+06 0.36 2.2256e+063.0509e+07 0.49 3.0509e+074.1823e+08 0.64 4.1823e+085.7333e+09 0.81 5.7333e+097.8594e+10 1 7.8594e+10步长为： 2.500000e-02四阶经典R-K法准确解误差0.33333 0.33333 00.013015 0.012246 0.000768940.040063 0.040015 4.82e-050.090037 0.09 3.6857e-050.16004 0.16 3.6723e-050.25004 0.25 3.6722e-050.36004 0.36 3.6722e-050.49004 0.49 3.6722e-050.64004 0.64 3.6722e-050.81004 0.81 3.6722e-051 1 3.6722e-05步长为： 1.000000e-02四阶经典R-K法准确解误差0.33333 0.33333 00.012256 0.012246 9.5673e-060.040016 0.040015 7.8252e-070.090001 0.09 6.6347e-070.16 0.16 6.6226e-070.25 0.25 6.6225e-070.36 0.36 6.6225e-070.49 0.49 6.6225e-070.64 0.64 6.6225e-070.81 0.81 6.6225e-071 1 6.6225e-07由结果可知，步长越小，结果越精确。
数值分析计算实习题列主元高斯消去法解线性方程组
数值分析计算实习题第5章解线性方程组的直接方法列主元高斯消去法解线性方程组。
书上的计算实习题1、2、3都要求用列主元高斯消去法解线性方程组，所以考虑写一个普适的程序来实现。
对于线性方程组Ax=b,程序允许用户从文件读入矩阵数据或直接在屏幕输入数据。
文件输入格式要求：(1)第一行为一个整数n (2<=n<=100),表示矩阵阶数。
(2)第2〜n+1行为矩阵A各行列的值。
(3 )第n+2〜n+n+2行为矩阵b各行的值。
屏幕输入：按提示输入各个数据。
输出：A、b、det(A)、列主元高斯消去计算过程、解向量X。
【算法说明】设有线性方程组Ax=b,其中设A为非奇异矩阵。
方程组的增广矩阵为[a,b] =a n…第1步(k=l)：首先在A的第一列中选取绝对值最大的元素％,作为第一步的主元素：如凜产如工0然后交换(A, b)的第1行与第1行元素，再进行消元计算。
设列主元素消去法已经完成第1步到第k・l步的按列选主元，交换两行，消元计算得到与原方程组等价的方程组A(k)x=b(k)用姐…哦盘…於？MB•■•■■[A,b]T[A（?b （*）] =■…盅•a"〉■唱”■■■••••...卅•第k步计算如下：对于k=l, 2, n-1 |時卜maxaf"(1)按列选主元：即确定t使(2)如果tHk,则交换［A, b］第t行与第k行元素。
(3)消元计算a ik J叫=-bgk + X・・Ji)a kk% <-%+叫ciy (i,j = k + l,…川)»〜乞+叫几 (j = k + l,…屮)消元乘数mik满足:M =吃 + 1,产•,)（4）回代求解（b厂工陶形）兀＜ ----- -- ---- ,（'=〃_1屮_2,・・、1）【程序】/*【普适列主元消去法解线性方程组】对于线性方程组：Ax=b输入：［选择屏幕直接输入］1・A的行阶数11(1 <=n <= 100)2.A的值3.b的值［选择读取文件］文件名(和主程序同级文件夹下)输出：1.A2.b3.det(A)4.解向量x#include <stdio.h>#include <stdlib.h>#include <math.h>double A[1O5][1O5],A_B[1O5][1O 习、b[105],x[105]; double detA;int njnaik = 1;〃读入数据void input(){int ij;char ch[20],name[100];FILE *f;printfC^—-\ii是否从文件读取数据(Y/N):”); scaufp%s：&ch);if(ch[O] = Y' || ch[0] = V){prmtfC'请输入文件名(包扌舌扩展名)：”)；scanfi(H%s,\name);f= fbpen(name；,i n);fscaiif(f/%d*\&n);fbr(i = 0;i < n;i ++)for(j = 0j<nj ++)fbr(i = 0;i < n;i ++)住canfg'%lf；&b[i]);pnntfC请输入A的阶数：”)；scanf(n%d %cT,&n);prmtfC请输入A的值：”)；fbr(i = 0;i < n;i ++)for(j = Oj<nj ++) scaiif(M%lf\&A[i][j]);printf("请输入b的值：”)；fbr(i = 0;i < n;i ++)}}〃计算行列式的值double det(double s[105][105],int m){int zj,k;double b[105][105],total = 0,r; /*b[N][N]ffl于存放,在矩阵s[N][N]中元素s[0]的余子式*/fbr(z = 0;z < m;z-H-){for(j = 0j<m-lj++)for(k = 0;k < m-l;k ++)if(k >= z)b[j][k] = s[j+l][k+l];elsebQ][k] = s[j+l][k];if(z % 2==0)r = s[0][z] * det(b4ii - 1); /*递归调用elser = (-1) * s[O]JXl * det(bjn - 1);total = total + r;}}else if(m = 2)total = s[0][0] * s[l][l] - s[0][l] * s[l][0];else if(m = 1)total = s[0][0];return total;}//输出A和b和det(A) void output_lQ{int ij;pnntffA = \ii n);for(i = 0;i < n;i ++){for。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
数值计算课程实习报告足球比赛中的吊门初射角问题——实习内容以及过程前言：吊射是球员在距离球门很远的一种射门方式，一般是看到对方门将站位靠前或者注意力有所分散时做吊射选择，由于足球在空中运行时间较长，可以给守门员反应的时间较长，需要精确的脚法才能成功，因此在职业比赛中出现吊射的情况很罕见，但吊射往往能起到出其不意，打破僵局的作用。
吊射的运动轨迹与弹道轨迹相似，可能的影响因素有球与球门的距离a，守门员与球门的距离b，球门高h，守门员最大接球高度H，球在空中飞行时间t，球出脚的初速度v，与水平方向的夹角alpha（初射角）。
问题的提出与分析考虑如下的因素：球与球门的距离为a，守门员与球门的距离为b，球门高h，守门员最大摸高H，球出脚的初速度为v，与水平方向的夹角为alpha（称为初射角）．给定，h=2.44m，H=3.20m，v=30m/s，重力加速度g=10m/s2提出问题：针对下列几组数据分别给出能吊门成功的相应初射角范围，画出最小及最大角度的运动轨迹。
要求精度在小数点后至少第3位。
(所给数据见计算结果)。
·a=6m，b=1m；·a=10m，b=3m；·a=20m，b=5m。
问题分析:1.先考虑最简单情形，即不考虑空气阻力等…，此时，球的运动轨迹是抛物线，如果守门员不动，总有合适的角度使吊门成功。
2.这不是求一个角度值，而是求一个范围！通常的思路是把问题整理成两个方程求根问题：一个方程是求吊门成功的最小角度，一个方程是求吊门成功的最大角度。
3.有可能落地弹入球门，要考虑反弹入门的情况。
直观分析:1.最简单情形，抛射体的运动轨迹为抛物线方程如下2.借助于使用方便的数学软件，可直观地看到各种初射角对应的抛射体运动的轨迹图形。
最简情形程序：v=30;g=10; h=2.44;H=3.2;a=6;b=1;l=a-b;L=a*1.01; %1.01表示进入门里 x=0:0.01:L; %对x 采样i1=floor(a/0.01)+1; i2=floor((a-b)/0.01)+1; %取整数 alpha=1.5368:0.00001:1.538; n=length(alpha);%弧度 for i=1:n;[y,tfinal]=paosheti1(x,alpha(i),v,g);tH=l/(v*cos(alpha(i))); %从射门到球到守门员位置的时间if y(i1)>=h | y(i2)<=H; success(i)=0; else success(i)=1; endp lot(l,H,'r+',a,h,'r+'), hold on,2222cos 1sin 2()tan 2cos x vt y vt gt g y x x xv αααα=⎧⎪⎨=-⎪⎩⇒=-plot(x,y), grid, hold offtitle(['足球比赛中的吊门','初射角=',num2str(alpha(i),6) ,...' 守门员的移动时间=',num2str(tH)]),pauseendfigure(2)plot(alpha,success)程序抛射体轨迹函数（算法）：function [y,t]=paosheti1(x,alpha,v,g)% 函数返回抛射体轨迹及第一第一次落地飞行时间y=x*tan(alpha)-x.^2*g/(2*v^2*(cos(alpha))^2);% 不考虑是直接进门还是落地后进门的轨迹t=2*v*sin(alpha)/g; % 第一次落地前最大飞行时间xmax=v*cos(alpha)*t; %第一次落地前最大飞行X向距离n=length(x);for i=1:nif y(i)<0 %判断是落地后进门xx=x(i)-xmax;y(i)=xx*tan(alpha)-xx.^2*g/(2*v^2*(cos(alpha))^2);%计算落地后进门的后半段轨迹endend运行结果:1.对于第一组数据，吊门成功的最小角度1.5369（为弧度，下同），对应的时间大约在4.9179秒，最大角度1.53791，对应的时间是5.0689秒；·a=6m，b=1m；最小角度：最大角度2.对于第二组数据，吊门成功的最小角度 1.51441，对应的时间大约在4.1403秒，最大角度1.5159，对应的时间大约在4.2526秒；·a=10m，b=3m；最小角度：最大角度：3.对于第三组数据，吊门成功的最小角度1.457，对应的时间大约在4.4033秒，最大角度1.4603，对应的时间大约是4.5343秒·a=20m，b=5m。
最小角度：最大角度：初步结果分析以及问题的再分析1.遵循由简单到一般的建模原则，先考虑简单情形，得到初步结果，以此为基础，发现问题、分析问题，找到求解思路，并逐步将问题一般化，甚至可以发现逐步一般化的顺序（还是由简到繁，先将哪些方面使之更一般化，如空气阻力、守门员移动）；2.结果有一定的合理性；3.从近似计算角度分析，在允许的精度范围内，如上的“作图——观察——调整”不失为一种求解方法；注：相比球与守门员及与球门的距离，注意守门员移动的时间，显然守门员有足够的时间移动，因此吊门是不会成功的！原因在于将问题假设得过于简单化、理想化了！修改假设应是下一步必须考虑的。
小结：最简情形之假设:1.不考虑空气阻力；2.不考虑守门员在球运行过程中的移动；3.球落地是完全弹性的，只考虑仅有一次触地反弹形成的吊门情况；4.只考虑越过守门员头顶的吊门，即出球点与守门员连成一线延伸到球门这样一个直线方向，不考虑从守门员侧面吊门的情况；5.将球看作是数学上的一个点；6.不考虑球的旋转，实际比赛时，旋转是很重要的！7.球的质量为一个单位。
有空气阻力的情形（一）——仅x 方向考虑空气阻力假设只考虑x 方向受空气阻力的影响；假设空气阻力与速度成正比，比例系数为k=0.4。
此时，x(t)满足如下的微分方程初值问题问题的解问题的解飞行轨迹：飞行到守门员位置的时间tH ：飞行到球门线位置的时间TH ：空气阻力的情形（一）（为了考虑时间不足，现在只针对第三组a=20;b=5时的情况，下同） 程序v=30;k=0.4;g=10;h=2.44;H=3.2;a=20;b=5;l=a-b;L=a*1.1;for alpha=1.2:0.001:1.3 %1.5425%pi/2-epsTh=-log(1-a*k/(v*cos(alpha)))/k;T=Th*1.2;t=0:0.01:T;1()cos (1)kt x t v e k α-=-21()sin 2y t vt gt α=-1()cos (1)1ln(1)cos ()kt x t v e kl t k k v x t lαα-⎧=-⎪⇒=--⎨⎪=⎩1()cos (1)1ln(1)cos ()kt x t v e ka t k k v x t a αα-⎧=-⎪⇒=--⎨⎪=⎩0x kx +=(0)0x =(0)cos x v α=[x,y]=paosheti2(t,alpha,v,k,g);TH=-log(1-l*k/(v*cos(alpha)))/k;plot(l,H,'r+',a,h,'r+'),hold on,plot(x,y),grid,hold offtitle(['足球比赛中的吊门','初射角=',num2str(alpha,6) ,...' 守门员的移动时间=',num2str(TH)]),pauseend空气阻力的情形（一）算法2：function [x,y]=paosheti2(t,alpha,v,k,g)x=v*cos(alpha)*(1-exp(-k*t))/k;y=v*sin(alpha)*t-g*t.^2/2;n=length(t);t0=2*v*sin(alpha)/g;%the time when the ball down to groundxt0=v*cos(alpha)*(1-exp(-k*t0))/k;vxt0=v*cos(alpha)*exp(-k*t0);vyt0=v*sin(alpha);vt0=sqrt(vxt0^2+vyt0^2);for i=1:nif t(i)>t0tt=t(i)-t0;x(i)=xt0+vt0*cos(alpha)*(1-exp(-k*tt))/k;y(i)=vt0*sin(alpha)*tt-g*tt^2/2;endendfunction [x,y]=paosheti2(t,alpha,v,k,g)%返回抛物体轨迹（考虑空气阻力）x=v*cos(alpha)*(1-exp(-k*t))/k;y=v*sin(alpha)*t-g*t.^2/2;n=length(t);t0=2*v*sin(alpha)/g; % 第一次落地前最大飞行时间xt0=v*cos(alpha)*(1-exp(-k*t0))/k; % 第一次落地前X向最大飞行距离vxt0=v*cos(alpha)*exp(-k*t0); % 第一次落地X向飞行速度vyt0=v*sin(alpha); % 第一次落地Y向飞行速度vt0=sqrt(vxt0^2+vyt0^2); % 第一次落地飞行速度大小for i=1:nif t(i)>t0 % 判断是落地后进门tt=t(i)-t0; % 计算落地后进门的后半段轨迹x(i)=xt0+vt0*cos(alpha)*(1-exp(-k*tt))/k;y(i)=vt0*sin(alpha)*tt-g*tt^2/2;endend空气阻力的情形（一）的结果及分析改进1前面结果有问题，反弹后的角度不应该是alpha 了，应该以落地时的情况计算出新反射角。
2.修改抛射体函数：将paosheti2(t,alpha,v,k,g)，换成paosheti22(t,alpha,v,k,g)。
tan()y x v alpha a v算法2修改如下得到算法22function [x,y]=paosheti22(t,alpha,v,k,g)x=v*cos(alpha)*(1-exp(-k*t))/k;y=v*sin(alpha)*t-g*t.^2/2;n=length(t);t0=2*v*sin(alpha)/g;%the time when the ball down to the groundxt0=v*cos(alpha)*(1-exp(-k*t0))/k;vxt0=v*cos(alpha)*exp(-k*t0);vyt0=v*sin(alpha);vt0=sqrt(vxt0^2+vyt0^2);alpha1=atan(vyt0/vxt0);for i=1:nif t(i)>t0tt=t(i)-t0;x(i)=xt0+vt0*cos(alpha1)*(1-exp(-k*tt))/k;y(i)=vt0*sin(alpha1)*tt-g*tt^2/2;endendfunction [x,y]=paosheti22(t,alpha,v,k,g)%返回抛物体轨迹（考虑空气阻力、反弹角度）x=v*cos(alpha)*(1-exp(-k*t))/k;y=v*sin(alpha)*t-g*t.^2/2;n=length(t);t0=2*v*sin(alpha)/g; % 第一次落地前最大飞行时间xt0=v*cos(alpha)*(1-exp(-k*t0))/k; % 第一次落地前X向最大飞行距离vxt0=v*cos(alpha)*exp(-k*t0); % 第一次落地X向飞行速度vyt0=v*sin(alpha); % 第一次落地Y向飞行速度vt0=sqrt(vxt0^2+vyt0^2); % 第一次落地飞行速度大小alpha1=atan(vyt0/vxt0); % 计算反弹角度for i=1:nif t(i)>t0 % 判断是落地后进门tt=t(i)-t0; % 计算落地后进门的后半段轨迹x(i)=xt0+vt0*cos(alpha1)*(1-exp(-k*tt))/k;y(i)=vt0*sin(alpha1)*tt-g*tt^2/2;endend空气阻力的情形（一）的进一步结果及分析改进针对第三组数据，计算的最小角度为1.2684，守门员移动时间为2.7836秒，最大角度是1.2807，时间是3.0032秒；最小角度：最大角度：结果仍有问题：反弹前后的两波高度一样；解决的办法是再考虑y 方向也有空气阻力。