西安交通大学计算方法B上机报告
西安交通大学汇编第二次上机实验报告范文
西安交通大学汇编第二次上机实验报告范文提交上机结果的模板文件第2次线上上机班级学号姓名1、循环程序设计-1（1）汇编、连接后的截图TODO：你的截图（必选）TODO：你的文字解释说明（可选）说明：mam某un得到某un.obj,某un.crf,某un.lt文件，通过link 某un得到某un.map文件，显示编译成功。
（2）.lt文件的截图，TODO：你的截图（必选）TODO：你的文字解释说明（可选）说明：通过mam对程序进行编译时生成.lt文件，通过notepad++打开.lt文件，并进行截图（3）反汇编的截图TODO：你的截图（必选）TODO：你的文字解释说明（可选）说明：在debug环境下执行u指令，显示出反汇编代码。
（4）在完成DS赋值后，立即显示各个寄存器的值TODO：你的截图（必选）TODO：你的文字解释说明（可选）说明：按单步t之后，下方出现MOVDS,A某，即下一条即将执行的指令为MOVDS,A某，再按一次t，此时DS被赋值，此时执行的r指令显示的就是DS赋值后各个寄存器的值。
（5）在进行计算前，显示数组M开始的n+2个字的内存值的截图（只能显示这n+2个字的内存值，多显示、少显示均扣分）TODO：你的截图（必选）TODO：你的文字解释说明（可选）说明：在debug环境下，执行d指令显示内存，由于不能显示其他字的内存值，所以只能一行一行截图，可以看到，此时内存值与程序初始定义值相同。
（6）执行完计算后，立即显示各个寄存器的值TODO：你的截图（必选）TODO：你的文字解释说明（可选）说明：先执行g指令到执行回到do系统指令，此时显然已经执行完运算，此时执行r指令就获得执行完运算后各个寄存器的值。
（7）执行完计算后，显示数组M开始的n+2个字的内存值的截图（只能显示这n+2个字的内存值，多显示、少显示均扣分）TODO：你的截图（必选）TODO：你的文字解释说明（可选）说明：执行d指令显示出内存值，由于要求是不能显示其他字的值，所以只能一行一行截图，可以看到此时内存值与期望结果相同。
计算方法与实习上机实验报告
计算方法与实习上机实验报告一、引言本文旨在介绍和展示我们在“计算方法”课程中的实习上机实验环节所完成的一些关键任务和所取得的成果。
该实验课程的目标是让我们更深入地理解和应用各种计算方法，并在实际操作中提高我们的编程和问题解决能力。
二、实验内容与目标实验的主要内容是利用各种计算方法解决实际数学问题。
我们被要求使用编程语言（如Python或Java）来实现和解决这些问题。
这些问题包括使用牛顿法求解平方根，使用蒙特卡洛方法计算圆周率，以及使用最优化方法求解函数的最小值等。
实验的目标不仅是让我们掌握计算方法的基本理论，更是要让我们能够在实际操作中运用这些方法。
我们需要在实习过程中，通过与同伴们合作，共同解决问题，提高我们的团队合作能力和问题解决能力。
三、实验过程与问题解决策略在实验过程中，我们遇到了许多问题，如编程错误、理解困难和时间压力等。
我们通过相互讨论、查阅资料和寻求教师帮助等方式，成功地解决了这些问题。
例如，在实现牛顿法求解平方根时，我们一开始对导数的计算和理解出现了一些错误。
但我们通过查阅相关资料和讨论，最终理解了导数的正确计算方法，并成功地实现了牛顿法。
四、实验结果与结论通过这次实习上机实验，我们不仅深入理解了计算方法的基本理论，还在实际操作中提高了我们的编程和问题解决能力。
我们的成果包括编写出了能有效求解平方根、计算圆周率和求解函数最小值的程序。
这次实习上机实验非常成功。
我们的团队不仅在理论学习和实践操作上取得了显著的进步，还在团队合作和问题解决方面积累了宝贵的经验。
这次实验使我们对计算方法有了更深的理解和认识，也提高了我们的编程技能和解决问题的能力。
五、反思与展望回顾这次实验，我们意识到在实验过程中，我们需要更好地管理我们的时间和压力。
在解决问题时，我们需要更有效地利用我们的知识和资源。
在未来，我们希望能够更加熟练地运用计算方法，并能够更有效地解决问题。
我们也希望能够将所学的计算方法应用到更广泛的领域中，如数据分析、科学研究和工业生产等。
西安交通大学算法上机实验报告
《计算机算法设计与分析》上机实验报告姓名：班级：学号：日期：2016年12月23日算法实现题3-14 最少费用购物问题★问题描述：商店中每种商品都有标价。
例如，一朵花的价格是2元，一个花瓶的价格是5元。
为了吸引顾客，商店提供了一组优惠商品价。
优惠商品是把一种或多种商品分成一组，并降价销售。
例如，3朵花的价格不是6元而是5元。
2个花瓶加1朵花的优惠价格是10元。
试设计一个算法，计算出某一顾客所购商品应付的最少费用。
★算法设计：对于给定欲购商品的价格和数量，以及优惠价格，计算所购商品应付的最少费用。
★数据输入：由文件input.txt提供欲购商品数据。
文件的第1行中有1个整数B（0≤B≤5），表示所购商品种类数。
在接下来的B行中，每行有3个数C,K和P。
C表示商品的编码（每种商品有唯一编码），1≤C≤999；K表示购买该种商品总数，1≤K≤5；P是该种商品的正常单价（每件商品的价格），1≤P≤999。
请注意，一次最多可购买5*5=25件商品。
由文件offer.txt提供优惠商品价数据。
文件的第1行中有1个整数S（0≤S≤99），表示共有S种优惠商品组合。
接下来的S行，每行的第1个数描述优惠商品组合中商品的种类数j。
接着是j个数字对（C，K），其中C是商品编码，1≤C≤999；K表示该种商品在此组合中的数量，1≤K≤5。
每行最后一个数字P （1≤P≤9999）表示此商品组合的优惠价。
★结果输出：将计算出的所购商品应付的最少费用输出到文件output.txt。
输入文件示例输出文件示例Input.txt offer.txt output.txt2 2 147 3 2 1 7 3 58 2 5 2 7 1 8 2 10解：设cost（a，b，c，d，e）表示购买商品组合（a，b，c，d，e）需要的最少费用。
A[k]，B[k]，C[k]，D[k]，E[k]表示第k种优惠方案的商品组合。
offer （m）是第m种优惠方案的价格。
西安交通大学数学建模上机实验报告
问题一某大型制药厂销售部门为了找出某种注射药品销量与价钱之间的关系,通过市场调查搜集了过去30个销售周期的销量及销售价钱的数据,如表.按照这些数据至少成立两个数学模型, 作出图形,比较误差。
问题分析：该问题是通过已知的过去30个销售周期的销量及销售价钱的 数据，来寻觅一个最能反映该药销量与价钱之间的函数曲 线。
在数学上归结为最佳曲线拟合问题。
大体思想：曲线拟合问题的提法:已知一组二维数据，即平面上的n 个点),x i i y （ i=1,2,3.....n ，i x 互不相同，寻求一个函数)(f y x =，使)(x f 在某中准则下与所有数据点最为接近，即曲线拟合得最好。
最小二乘法是解决曲线拟合最常常利用的方式.大体思路：1122 ()()()()m m f x a r x a r x a r x =+++令其中rk(x) 是事前选定的一组函数，ak 是待定系数(k=1,2,…,m,m ＜n), 拟合准则是使n 个点(xi,yi) (i=1,2…,n)，与y=f(xi)的距离 的平方和最小，称最小二乘法准则。
一、系数的肯定22111 (,,)[()]n nm ii i i i J a a f x y δ====-∑∑记求m a a ,,1 使得使J 达到最小.0 (1,,)kJ k m a ∂==∂ 取得关于 m a a ,,1 的线性方程组：11111()[()]0 ()[()]0nmi k k i i i k n mm i k k i i i k r x a r x y r x a r x y ====⎧-=⎪⎪⎪⎨⎪⎪-=⎪⎩∑∑∑∑ 1 ,,().m a a f x 解出，即得散点图： 程序： x=[,,,,,,,,,,,,,,,,,,,,,,,,,,,,,]; y=[,,,,,,,,,,,,,,,,,,,,,,,,,,,,,]; plot(x,y,'r.')通过观察，结合实际情形。
西安交大计算方法上机报告
从而得到计算的公式：
j 1, 2,..., n 1 j 1 j l i1 i 2,3,..., n i1 11 i 1 lik ukj j i , i 1,..., n, i 2,3,.., n ij ij k 1 i 1 1 l ( lkt ti ) k i 1,..., n, i 2,3,.., n ki ki t 1 ii
计算方法 上机实习题目报告
班级：材料 s3076 班 姓名：丁明帅 学号：3113305029
1：计算 S
100000
k 1
1 ,要求误差小于 106 ,给出实现算法。 k2
【实现思路】
设当 k 值为 i 时 S-Si<10-6,则只需要
S Si
100000 1
k i
1
2
1 l32 ln 2 1 ln 3 1
4 / 46
12 22
1n 2 n
nn
因此有
1 j 2 j 0 jj 0 0
ij li1
li 2
li ,i -1
根据定义知道 L 矩阵的斜对角线上的值都为 1，且 L 矩阵的第一行与原矩阵 A 的第一行 相同，因此可以根据公式先计算 U 的第一行，然后计算 L 的第一列；以后的第 i 步先计算 U 的第 i 行， 然后计算 L 的第 i 列 （U 的第 n 行不作计算） 。 然后把最后的 U 的第 n 行计算出来。
【算法依据】
列主元高斯消元法的过程可以将方程组系数简化为系数矩阵与 b 矩阵， 从而利用方程组 对系数扩展矩阵进行消元。 在消元的过程中矩阵的行向量之间可以变换， 但列向量不能变化。 在进行压缩矩阵的求解中还需要人为的将因调整行向量所导致的列向量的变化调整回来。 Matlab 对于矩阵的处理非常容易，结合 for 循环语句可以实验本题大规模方程组的求解。
西安交大 计算方法B上机作业
计算方法（B ）上机作业一、三次样条拟合某通信公司在一次施工中，需要在水面宽度为20米的河沟底部沿直线走向铺设一条沟底光缆。
在铺设光缆之前需要对沟底的地形进行初步探测，从而估计所需光缆的长度，为工程预算提供依据。
已探测到一组等分点位置的深度数据(单位:米)如下表所示：(1)(2)估算所需光缆长度的近似值，并作出铺设河底光缆的曲线图； 解：1、算法实现的思想及依据题目（1）为曲线拟合问题多项式插值、分段插值和最小二乘法。
多项式插值，随着插值数据点的数目增多，误差也会随之增大，因此不选用。
最小二乘法适于数据点较多的场合，在此也不适用。
故选用分段插值。
分段插值又分为分段线性插值、分段二次插值、三次样条插值及更高阶的多项式插值。
由本题的物理背景知，光缆正常工作时各点应该是平滑过渡，因此至少选用三次样条插值法。
对于更高阶的多项式插值，由于“龙格现象”而不选用。
题目（2）求光缆长度，即求拟合曲线在0到20的长度，对弧长进行积分即可。
光缆长度的第一型线积分表达式为190k kk l +==∑⎰。
2、算法实现的结构参考教材给出的SPLINEM 算法和TTS 算法，在选定边界条件和选定插值点等距分布后，可以先将数据点的二阶差商求出来并赋值给右端向量d ，再根据TSS 解法求解M 。
光缆长度的第一型线积分表达式为190k kk l +==∑⎰。
3、程序运行结果及分析图1.1三种边界条件下三次样条插值图1.2光缆长度4、MATLAB代码：1）自己编程实现代码clear;clc;I=input('你想使用第几种边界条件？请输入1、2、3之一: ');x=0:20;y=[9.01 8.96 7.96 7.97 8.02 9.05 10.13 11.18 12.26 13.28 13.32 12.61 11.29 10.22 9.15 7.90 7.95 8.86 9.81 10.8 10.93];plot(x,-y,'k.','markersize',15)%y为深度，取负号hold on%% 计算一阶差商y1=ones(1,21);for i=2:1:21y1(i)=(y(i)-y(i-1))/(x(i)-x(i-1));end%% 计算二阶差商y2=ones(1,21);for i=3:1:21y2(i)=(y1(i)-y1(i-1))/(x(i)-x(i-2));end%% 计算三阶差商y3=ones(1,21);for i=4:1:21y3(i)=(y2(i)-y2(i-1))/(x(i)-x(i-3));end%% 选择边界条件（I）if I==1d(1)=0;d(21)=0;a(21)=0;c(1)=0;% 第一个点和最后一个点的二阶差商为0 endif I==2d(1)=6*y1(1);d(21)=-6*y1(21);a(1)=1;c(1)=1;endif I==3d(1)=-12*y3(1);d(21)=-12*y3(21);a(21)=-2;c(1)=-2;%endfor i=2:20d(i)=6*y2(i+1);end%% 构造带状矩阵求解（追赶法）b=2*ones(1,21);a=0.5*ones(1,21);%a(21)=-2;c=0.5*ones(1,21);%c(1)=-2;u(1)=b(1);r(1)=c(1);%% 追yz(1)=d(1);for i=2:21l(i)=a(i)/u(i-1);u(i)=b(i)-l(i)*r(i-1);r(i)=c(i);yz(i)=d(i)-l(i)*yz(i-1);end%% 赶xg(21)=yz(21)/u(21);for i=20:-1:1xg(i)=(yz(i)-r(i)*xg(i+1))/u(i);endM=xg;%%所有点的二阶导数值%% 求函数表达式并积分t=1;h=1;N=1000x1=0:20/(N-1):20;length=0;for i=1:Nfor j=2:20if x1(i)<=x(j)t=j;break;elset=j+1;endendf1=x(t)-x1(i);f2=x1(i)-x(t-1);S(i)=(M(t-1)*f1^3/6/h+M(t)*f2^3/6/h+(y(t-1)-M(t-1)*h^2/6)*f1+(y(t)-M(t)*h^2/6)* f2)/h;Sp(i)=-M(t-1)*f1^2/2/h+M(t)*f2^2/2/h+(y(t)-y(t-1))/h-(M(t)-M(t-1))*h/6;length(i+1)=sqrt(1+Sp(i)^2)*(20/(N-1))+length(i);%第一类线积分endfigure(1);plot(x1,-S,'r-')%深度曲线griddisp(['第',num2str(I),'种边界条件下长度',num2str(length(N+1)),'米'])axis fill;xlabel('测点/米');ylabel('深度/米');title('三次样条曲线拟合');legend('数据点','拟合曲线',3);二、最小二乘近似假定某天的气温变化记录如下表所示，试用数据拟合的方法找出这一天的气温变化的规律;试计算这一天的平均气温,并试估计误差。
西安交大计算方法B2017大作业
计算方法B上机报告姓名：学号：班级：学院：任课教师：2017年12月29日题目一：1.1题目内容某通信公司在一次施工中，需要在水面宽度为20米的河沟底部沿直线走向铺设一条沟底光缆。
在铺设光缆之前需要对沟底的地形进行初步探测，从而估计所需光缆的长度，为工程预算提供依据。
已探测到一组等分点位置的深度数据(单位:(1)(2)估算所需光缆长度的近似值，并作出铺设河底光缆的曲线图；1.2 实现题目的思想及算法依据首先在题目（1）中要实现的是数据的拟合，显然用到的是我们在第三章中数据近似的知识内容。
多项式插值时，这里有21个数据点，则是一个20次的多项式，但是多项式插值随着数据点的增多，会导致误差也会随之增大，插值结果会出现龙格现象，所以不适用于该题目中点数较多的情况。
为了避免结果出现大的误差，同时又希望尽可能多地使用所提供的数据点，提高数据点的有效使用率，这里选择分段插值方法进行数据拟合。
分段插值又可分为分段线性插值、分段二次插值和三次样条插值。
由于题目中所求光缆的现实意义，而前两者在节点处的光滑性较差，因此在这里选择使用三次样条插值。
根据课本SPLINEM 算法和TSS 算法，采用第三种真正的自然边界条件，在选定边界条件和选定插值点等距分布后，可以先将数据点的二阶差商求出并赋值给右端向量d ，再根据TSS 解法求解三对角线线性方程组从而解得M 值。
求出M 后，对区间进行加密，计算200个点以便于绘图以及光缆长度计算。
对于问题（2），使用以下的公式:20=()L f x ds ⎰20(f x =⎰191(k kk f x +==∑⎰1.3 算法结构1. For n i ,,2,1,0⋅⋅⋅=1.1 i i M y ⇒2. For 2,1=k2.1 For k n n i ,,1, -=2.1.1 i k i i i i M x x M M ⇒----)/()(13. 101h x x ⇒-4. For 1-,,2,1n i =4.1 11++⇒-i i i h x x4.2 b a c c h h h i i i i i i ⇒⇒-⇒+++2;1;)/(11 4.3 i i d M ⇒+165. 0000;;c M d M d n n ⇒⇒⇒λn n n b a b ⇒⇒⇒2;;20μ6. 1111,γμ⇒⇒d b7. For m k ,,3,2 = ! 获取M 的矩阵元素个数，存入m7.1 k k k l a ⇒-1/μ 7.2 k k k k c l b μ⇒⋅-1- 7.3 k k k k l d γγ⇒⋅-1- 8. m m m M ⇒μγ/9. For 1,,2,1 --=m m k9.1 k k k k k M M c ⇒⋅-+μγ/)(110. k ⇒1 ! 获取x 的元素个数存入s 11. For 1,,2,1-=s i11.1 if i x x ≤~then k i ⇒；break else k i ⇒+112. xx x x x x h x x k k k k ˆ~;~;11⇒-⇒-⇒--- y h x h M y x h M y x M x M k k k k k k ~/]ˆ)6()6(6ˆ6[2211331⇒-+-++---1.4 matlab 源程序n=20; x=0:n;y=[9.01 8.96 7.96 7.97 8.02 9.05 10.13 11.18 12.26 13.28 13.32 12.61 11.29 10.22 9.15 7.90 7.95 8.86 9.81 10.80 10.93];M=y; %用于存放差商，此时为零阶差商 h=zeros(1,n+1); c=zeros(1,n+1); d=zeros(1,n+1); a=zeros(1,n+1); b=2*ones(1,n+1); h(2)=x(2)-x(1);for i=2:n %书本110页算法SPLINEM h(i+1)=x(i+1)-x(i);c(i)=h(i+1)/(h(i)+h(i+1)); a(i)=1-c(i); enda(n+1)=-2; %计算边界条件c(0),a(n+1),采用的是第三类边界条件 c(1)=-2;for k=1:3 %计算k 阶差商 for i=n+1:-1:k+1M(i)=(M(i)-M(i-1))/(x(i)-x(i-k)); endif(k==2) %计算2阶差商 d(2:n)=6*M(3:n+1); %给d 赋值 endif(k==3)d(1)=(-12)*h(2)*M(4); %计算边界条件d(0),d(n),采用的是第三类边界条件 d(n+1)=12*h(n+1)*M(n+1); end endl=zeros(1,n+1); r=zeros(1,n+1); u=zeros(1,n+1); q=zeros(1,n+1); u(1)=b(1); r(1)=c(1); q(1)=d(1);for k=2:n+1 %利用书本49页算法TSS求解三对角线性方程组r(k)=c(k);l(k)=a(k)/u(k-1);u(k)=b(k)-l(k)*r(k-1);q(k)=d(k)-l(k)*q(k-1);endp(n+1)=q(n+1)/u(n+1);for k=n:-1:1p(k)=(q(k)-r(k)*p(k+1))/u(k);endfprintf('三对角线性方程组的解为：');disp(p);%求拟合曲线x1=0:0.1:20; %首先对区间进行加密,增加插值点n1=10*n;x2=zeros(1,n1+1);x3=zeros(1,n1+1);s=zeros(1,n1+1);for i=1:n1+1for j=1:nif x1(i)>=x(j)&&x1(i)<=x(j+1) %利用书本111页算法EVASPLINE求解拟合曲线s(x)h(j+1)=x(j+1)-x(j);x2(i)=x(j+1)-x1(i);x3(i)=x1(i)-x(j);s(i)=(p(j).*(x2(i)).^3/6+p(j+1).*(x3(i)).^3/6+(y(j)-p(j).*((h(j+1)).^2/6)).*x2( i)+...(y(j+1)-p(j+1).*(h(j+1)).^2/6).*x3(i))/h(j+1);endendendplot(x,-y,'x') %画出插值点hold onplot(x1,-s) %画出三次样条插值拟合曲线hold ontitle('三次样条插值法拟合电缆曲线');xlabel('河流宽度/m');ylabel('河流深度/m');Length=0;for i=1:n1L=sqrt((x1(i+1)-x1(i))^2+(s(i+1)-s(i))^2); %计算电缆长度Length=Length+L;endfprintf('电缆长度(m)=');disp(Length);1.5 结果与说明铺设海底光缆的曲线如图1.1所示图1. 1三次样条插值法拟合海底光缆曲线由上图可以看出，所得到的曲线光滑，能够较好得反映实际的河沟底部地势形貌。
西安交通大学计算方法B完整版
第一章绪论1.1数值计算现代科学的发展，已导致科学与技术的研究从定性前进到定量，尤其是现代数字计算机的出现及迅速发展，为复杂数学问题的定量研究与解决，提供了强有力的基础。
通常我们面对的理论与技术问题，绝大多数都可以从其物理模型中抽象出数学模型，因此，求解这些数学模型已成为我们面临的重要任务。
一、本课程的任务：寻求解决各种数学问题的数值方法——如何将高等数学的问题回归到初等数学（算术）的方法求解——了解计算的基础方法，基本结构（否则只须知道数值软件）——并研究其性质。
立足点：面向数学——解决数学问题面向计算机——利用计算机作为工具充分发挥计算机的功能，设计算法，解决数学问题例如：迭代法、并行算法二、问题的类型1、离散问题：例如，求解线性方程组bAx=——从离散数据：矩阵A和向量b，求解离散数据x；2、连续问题的离散化处理：例如，数值积分、数值微分、微分方程数值解；3、离散问题的连续化处理：例如，数据近似，统计分析计算；1.2数值方法的分析在本章中我们不具体讨论算法，首先讨论算法分析的基础——误差。
一般来讲，误差主要有两类、三种（对科学计算）：1）公式误差——“截断误差”，数学↔计算，算法形成——主观（人为）：数学问题-数值方法的转换，用离散公式近似连续的数学函数进行计算时，一般都会发生误差，通常称之为“截断误差”；——以后讨论2）舍入误差及输出入误差——计算机，算法执行——客观（机器）：由于计算机的存储器、运算器的字长有限，在运算和存储中必然会发生最末若干位数字的舍入，形成舍入误差；在人机数据交换过程中，十进制数和二进制数的转换也会导致误差发生，这就是输入误差。
这两种误差主要是由于计算机的字长有限，采用浮点数系所致。
首先介绍浮点数系一、计算机上的运算——浮点运算面向计算机设计的算法，则先要讨论在计算机上数的表示。
科学记数法——浮点数：约定尾数中小数点之前的数全为零，小数点后第一个数不能为零。
目前，一般计算机都采用浮点数系，一个存储单元分成首数和尾数：首数l 尾数(位)其中首数存放数的指数（或“阶”）部分，尾数存放有效数字。
西安交通大学传热学上机报告材料-墙角导热数值分析报告
实用文档传热大作业二维导热物体温度场的数值模拟：璇班级：能动A02学号：10031096一．物理问题有一个用砖砌成的长方形截面的冷空气通道，其截面尺寸如下图所示，假设在垂直于纸面方向上用冷空气及砖墙的温度变化很小，可以近似地予以忽略。
在下列两种情况下试计算：（1）砖墙横截面上的温度分布；（2）垂直于纸面方向的每米长度上通过砖墙的导热量。
第一种情况：外壁分别均与地维持在0℃及30℃；第二种情况：外壁均为第三类边界条件，且已知：t ∞1=30℃,ℎ1=10wm2∙℃t ∞2=10℃,ℎ2=4wm2∙℃砖墙的导热系数λ=0.53 Wm∙℃二．数学描写由对称的界面必是绝热面，可取左上方的四分之一墙角为研究对象，该问题为二维、稳态、无热源的导热问题，其控制方程和边界条件如下：ðt2ðx2+ðt2ðy2=0边界条件（情况一） t(x,0)=30 0≤x≤1.5t(0,y)=30 0≤y≤1.1t(0.5,y)=0 0.5≤y≤1.1t(x,0.5)=0 0.5≤x≤1.5ðt(1.5,y)=0 0≤y≤0.5ðy∂t(x,1.1)=0 0≤x≤0.5ℎ(t−t f1) x=0,0≤y≤1.11=ℎ(t−t f2) x=0.5,0.5≤y≤1.12ℎ(t−t f1) y=0,0≤x≤1.51ℎ(t−t f2) y=0,0.5≤x≤21.50 0≤y≤0.5=0 0≤x≤0.5∂x三．网格划分网格划分与传热学实验指导书中“二维导热物体温度场的电模拟实验”一致，如下图所示：四．方程离散对于节点，离散方程t[i][j]=0.25*(t[i+1][j]+t[i-1][j]+t[i][j+1]+t[i][j-1])对于边界节点，则应对一、二两种情况分开讨论：情况一:绝热平直边界点：t[15][j]=0.25*(2*t[14][j]+t[15][j-1]+t[15][j+1])1≤j≤4t[i][11]=0.25*(2*t[i][10]+t[i-1][11]+t[i+1][11]) 1≤i≤4外等温边界点：t[i][j]=30等温边界点：t[i][j]=0情况二：（Bi1,Bi2为网格Bi数，Bi1=ℎ1∆xλ Bi2=ℎ2∆xλ）绝热平直边界点：t[15][j]=0.25*(2*t[14][j]+t[15][j-1]+t[15][j+1])1≤j≤4t[i][11]=0.25*(2*t[i][10]+t[i-1][11]+t[i+1][11]) 1≤i≤4外侧对流平直边界：t[i][0]=(2*t[i][1]+t[i+1][0]+t[i-1][0]+2*Bi1*tf1)/(2*Bi1+4) 1≤i≤14t[0][j]=(2*t[1][j]+t[0][j+1]+t[0][j-1]+2*Bi1*tf1)/(2*Bi1+4) 1≤j≤10侧对流平直边界：t[i][5]=(2*t[i][4]+t[i+1][5]+t[i-1][5]+2*Bi2*tf2)/(2*Bi2+4) 6≤i≤14t[5][j]=(2*t[4][j]+t[5][j+1]+t[5][j-1]+2*Bi2*tf2)/(2*Bi2+4) 6≤j≤10特殊点：a点t[15][0]=(t[14][0]+t[15][1]+tf1*Bi1)/(Bi1+2)b点t[15][5]=(t[14][5]+t[15][4]+tf2*Bi2)/(Bi2+2)c点t[5][5]=(2*t[4][5]+2*t[5][4]+t[5][6]+t[6][5]+3*Bi2*tf2)/(2*Bi2+6) d点t[5][11]=(t[5][10]+t[4][11]+tf2*Bi2)/(Bi2+2)e点t[0][11]=(t[0][10]+t[1][11]+tf1*Bi1)/(Bi1+2)f点t[0][0]=(t[0][1]+t[1][0]+tf1*Bi1*2)/(2*Bi1+2)五．编程思路及流程图编程思路为设定两个二维数组t[i][j]、ta[i][j]分别表示本次迭代和上次迭代各节点的温度值，iter表示迭代进行的次数,daore_in、daore_out分别表示外边界的散热量。
西安交通大学 计算方法上机报告
计算方法上机报告姓名：学号：班级：机械硕4002上课班级：02班说明：本次上机实验使用的编程语言是Matlab 语言，编译环境为MATLAB 7.11.0，运行平台为Windows 7。
1. 对以下和式计算：∑∞⎪⎭⎫ ⎝⎛+-+-+-+=0681581482184161n n n n S n，要求：① 若只需保留11个有效数字，该如何进行计算； ② 若要保留30个有效数字，则又将如何进行计算；（1） 算法思想1、根据精度要求估计所加的项数，可以使用后验误差估计，通项为： 142111416818485861681n n n a n n n n n ε⎛⎫=---<< ⎪+++++⎝⎭； 2、为了保证计算结果的准确性，写程序时，从后向前计算； 3、使用Matlab 时，可以使用以下函数控制位数： digits(位数)或vpa(变量，精度为数)（2）算法结构1. ;0=s⎪⎭⎫⎝⎛+-+-+-+=681581482184161n n n n t n； 2. for 0,1,2,,n i =⋅⋅⋅ if 10m t -≤end;3. for ,1,2,,0n i i i =--⋅⋅⋅;s s t =+（3）Matlab源程序clear; %清除工作空间变量clc; %清除命令窗口命令m=input('请输入有效数字的位数m='); %输入有效数字的位数s=0;for n=0:50t=(1/16^n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6));if t<=10^(-m) %判断通项与精度的关系break;endend;fprintf('需要将n值加到n=%d\n',n-1); %需要将n值加到的数值for i=n-1:-1:0t=(1/16^i)*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6));s=s+t; %求和运算ends=vpa(s,m) %控制s的精度（4）结果与分析当保留11位有效数字时，需要将n值加到n=7，s =3.1415926536；当保留30位有效数字时，需要将n值加到n=22,s =3.14159265358979323846264338328。
数值计算方法上机实验报告
数值计算方法上机实验报告实验目的：复习和巩固数值计算方法的基本数学模型，全面掌握运用计算机进行数值计算的具体过程及相关问题。
利用计算机语言独立编写、调试数值计算方法程序，培养学生利用计算机和所学理论知识分析解决实际问题的能力。
上机练习任务：利用计算机基本C 语言编写并调试一系列数值方法计算通用程序，并能正确计算给定题目，掌握调试技能。
掌握文件使用编程技能，如文件的各类操作，数据格式设计、通用程序运行过程中文件输入输出运行方式设计等。
一、各算法的算法原理及计算机程序框图1. 列主元高斯消去法算法原理：高斯消去法是利用现行方程组初等变换中的一种变换，即用一个不为零的数乘一个方程后加只另一个方程，使方程组变成同解的上三角方程组，然后再自下而上对上三角方程组求解。
列选住院是当高斯消元到第k 步时，从k 列的kk a 以下（包括kk a ）的各元素中选出绝对值最大的，然后通过行交换将其交换到kk a 的位置上。
交换系数矩阵中的两行（包括常数项），只相当于两个方程的位置交换了，因此，列选主元不影响求解的结果。
●源程序：#define N 200#include "stdio.h"#include "math.h"FILE *fp1,*fp2;void LZ(){int n,i,j,k=0,l;double d,t,t1;static double x[N],a[N][N];fp1=fopen("a1.txt","r");fp2=fopen("b1.txt","w");fscanf(fp1,"%d",&n);for(i=0;i<n;++i)for(j=0;j<=n;++j){fscanf(fp1,"%lf",&a[i][j]);}{d=a[k][k];l=k;i=k+1;do{if(fabs(a[i][k])>fabs(d)) /*选主元*/{d=a[i][k];l=i;}i++;}while(i<n);if(d==0){printf("\n输入矩阵有误！\n");}else{ /*换行*/if(l!=k){for(j=k;j<=n;j++){t=a[l][j];a[l][j]=a[k][j];a[k][j]=t;}}}for(j=k+1;j<=n;j++) /*正消*/ a[k][j]/=a[k][k];for(i=k+1;i<n;i++)for(j=k+1;j<=n;j++)a[i][j]-=a[i][k]*a[k][j];k++;}while(k<n);if(k!=0){for(i=n-1;i>=0;i--) /*回代*/ {t1=0;for(j=i+1;j<n;j++)t1+=a[i][j]*x[j];x[i]=a[i][n]-t1;}for(i=0;i<n;i++)fprintf(fp2,"\n 方程组的根为x[%d]=%lf",i+1,x[i]); fclose(fp1); fclose(fp2); }main() { LZ(); }● 具体算例及求解结果：用列选主元法求解下列线性方程组⎪⎩⎪⎨⎧=++=++=-+28x x 23x 2232832321321321x x x x x x 输入3 输出结果：方程组的根为x[1]=6.0000001 2 -3 8 方程组的根为x[2]=4.000000 2 1 3 22 方程组的根为x[3]=2.000000 3 2 1 28● 输入变量、输出变量说明：输入变量：ij a 系数矩阵元素，i b 常向量元素 输出变量：12,,n b b b 解向量元素2. 杜里特尔分解法解线性方程● 算法原理：求解线性方程组Ax b =时，当对A 进行杜里特尔分解，则等价于求解LUx b =，这时可归结为利用递推计算相继求解两个三角形（系数矩阵为三角矩阵）方程组，用顺代，由Ly b =求出y ，再利用回带，由Ux y =求出x 。
西安交大计算方法
西安交通大学计算方法上机实验班级：(xxx)姓名：(xxx)学号：21116010041.按两种顺序计算y，哪个接近真值？Y = 1000 + + + … +用java 语言编写：public class Add {public static void main(String[] args){double s=0,y=1000;for(double a=1001.0;a<=2000.0;a++){y+=1.0/a;}for(double a=2000.0;a>=1001.0;a--){s+=1.0/a;}s=s+1000;System.out.println("正序和"+s);System.out.println("逆序和"+y);}}运行结果：结论：显然假设是double类型的数据时，先算大数的过程吃掉了末尾的小数被进位所埋没，导致了大数吃小数的误差，按从小到大（从右向左）的计算顺序所得的结果与真值相近，而按从大到小（从左到右）的计算顺序所得的结果与真值的误差较大。
1-18.设(x) = 1 + x + + + … + , 计算(-5)和1/(5)，哪个接近？解法一：用JAVA 语言编写：public class second{ public static void main(String[] args){double s1=1 ,s2=1;double e=1,sum=1; //e的初值为1，sum用来存放n！int a=1;while(sum<Math.pow(10, 1000000)){sum=a*sum;e=1.0/sum+e;a++;}double b=1.0/(e*e*e*e*e);System.out.println("较为精确的值1/e^5="+b);for(int i=1;i<=24;i++){s1+=cimi1(i);s2+=cimi2(i);}s1=1.0/s1;System.out.println("1/S24(5)="+s1);System.out.println("S24(-5)="+s2);}public static double cimi1(int ai){double xi=1;for(int i=ai;i>=1;i--){xi=xi*(5.0/i);}return xi;}public static double cimi2(int ai){double xi=1;for(int i=ai;i>=1;i--){xi=xi*(-5.0/i);}return xi;}}运行结果：解法二：用matlab编程并运行，如下：(1)计算(-5)运行结果如下：(2)计算1/(5)运行结果如下：而的真是结果为0.006737946比较得1/(5)的计算结果与真实值更接近解法三：也可以用C++编写：#include "stdafx.h"#include"stdio.h"#include "iostream"using namespace std;int main(int argc, char* argv[]){ int func1(int );double func2(int);double y=0;int i;for(i=1;i<25;i++){ int z=func1(i);double e=func2(i);y+=z/e;}cout<<"----------------------------------------"<<endl;cout<<"1/S(5)的运算结果是："<<" "<<1.0/(y+1)<<endl;cout<<"----------------------------------------"<<endl;return 0;}int func1(int x){int y=1;int k;for (k=0;k<x;k++)y*=5;return y;}double func2(int n){double y=1;int j;for (j=1;j<=n;j++)y*=j;return y;}运行结果如下图：结论：通过比较上述的几种编程结果，可以看出1/S(5),更接近真实值，而且用matlab更为简便，可以直接利用函数库，并可以轻松的嵌入秦九韶算法，大大减少运算量和时间。
计算方法b课程设计
计算方法b课程设计一、课程目标知识目标：1. 学生能理解并掌握计算方法B中的基本概念，如插值法、数值积分和微分方程的数值解法。
2. 学生能够运用适当的计算方法解决实际问题，并对结果进行分析和评价。
3. 学生能够掌握计算方法B中涉及的数学原理和算法步骤，与高年级数学课程内容相衔接。
技能目标：1. 学生能够运用计算软件（如MATLAB）实现计算方法B中的算法，解决相关数学问题。
2. 学生通过小组合作，培养解决复杂问题的能力，学会在实际情境中运用计算方法B。
3. 学生能够通过实际案例，分析数据，撰写报告，提高计算思维和解决问题的技能。
情感态度价值观目标：1. 学生对计算方法B产生兴趣，认识到其在工程、自然科学等领域的重要应用。
2. 学生在学习过程中，培养严谨的科学态度，树立正确的数学观念。
3. 学生通过小组合作，培养团队协作精神，增强沟通与交流能力。
课程性质：本课程为专业基础课程，旨在帮助学生建立扎实的计算方法基础，提高解决实际问题的能力。
学生特点：学生处于大学二年级，具备一定的数学基础和编程能力，但对计算方法B的了解有限。
教学要求：结合学生特点，注重理论与实践相结合，强化上机操作和实际案例分析，提高学生的实际应用能力。
在教学过程中，关注学生的个别差异，因材施教，确保每位学生都能达到课程目标。
通过课程学习，使学生在知识、技能和情感态度价值观方面取得具体的学习成果，为后续课程和实际工作打下坚实基础。
二、教学内容本课程教学内容主要包括以下几部分：1. 插值法：介绍拉格朗日插值、牛顿插值、分段插值等方法，分析各种插值法的优缺点及适用范围。
2. 数值积分：讲解梯形法、辛普森法、高斯求积等数值积分方法，探讨数值积分的误差分析。
3. 微分方程的数值解法：介绍欧拉法、改进欧拉法、龙格-库塔法等求解常微分方程初值问题的数值方法。
4. 线性方程组的数值解法：讲解高斯消元法、LU分解、迭代法等方法，分析算法的稳定性和收敛性。
计算方法实验上机报告(完整版)
简单迭代法#include<stdio.h>#include<math.h>#define x0 3.0#define MAXREPT 1000#define EPS 1E-6#define G(x) pow(12*x+sin(x)-1,1.0/3)void main(){int i;double x_k=x0,x_k1=x0;printf("k\txk\n");for(i=0;i<MAXREPT;i++){printf("%d\t%g\n",i,x_k1);x_k1=G(x_k);if (fabs(x_k1-x_k)<EPS){printf("THE ROOT IS x=%g,k=%d\n",x_k1,i);return;}x_k=x_k1;}printf("AFTER %d repeate,no solved.\n",MAXREPT);}结果牛顿迭代法一#include<stdio.h>#include<math.h>#define x0 3.0#define MAXREPT 1000#define EPS 1E-6#define G(x) x-(pow(x,3)-sin(x)-12*x+1)/(3*pow(x,2)-cos(x)-12) void main(){int i;double x_k=x0,x_k1=x0;printf("k\txk\n");for(i=0;i<MAXREPT;i++){printf("%d\t%g\n",i,x_k1);x_k1=G(x_k);if (fabs(x_k1-x_k)<EPS){printf("THE ROOT IS x=%g,k=%d\n",x_k1,i);return;}x_k=x_k1;}printf("AFTER %d repeate,no solved.\n",MAXREPT);}结果埃特金加速法#include<stdio.h>#include<math.h>#define x0 3.0#define MAXREPT 1000#define EPS 1E-6#define G(x) (pow(x,3)-sin(x)+1)/12void main(){int i;double x1=x0,x2=x0;double z,y;printf("k\tx1\tx2\txk\n");for(i=0;i<MAXREPT;i++){if(i==0)printf("%d\t\t\t%g\n",i,x2);elseprintf("%d\t%g\t%g\t%g\n",i,y,z,x2);y=G(x1);z=G(y);x2=z-((z-y)*(z-y))/(z-2*y+x1);if (fabs(x2-x1)<EPS){printf("THE ROOT IS x=%g,k=%d\n",x2,i);return;}x1=x2;}printf("AFTER %d repeate,no solved.\n",MAXREPT);} 结果牛顿迭代法二#include<stdio.h>#include<math.h>#define x0 1.5#define MAXREPT 1000#define EPS 1E-6#define G(x) x-(pow(x,3)+pow(x,2)-3*x-3)/(3*pow(x,2)+2*x-3) void main(){int i;double x_k=x0,x_k1=x0;printf("k\txk\n");for(i=0;i<MAXREPT;i++){printf("%d\t%g\n",i,x_k1);x_k1=G(x_k);if (fabs(x_k1-x_k)<EPS){printf("THE ROOT IS x=%g,k=%d\n",x_k1,i);return;}x_k=x_k1;}printf("AFTER %d repeate,no solved.\n",MAXREPT);}结果弦截法#include<stdio.h>#include<math.h>#define x0 0#define x1 1.5#define MAXREPT 1000#define EPS 1E-6#define G(x) pow(x,3)+pow(x,2)-3*x-3void main(){int i;double x_k=x0,x_k1=x1,x_k2=0;double y,z;printf("k\txk\n");for(i=0;i<MAXREPT;i++){printf("%d\t%g\n",i,x_k2);y=G(x_k);z=G(x_k1);x_k2=x_k1-(z*(x_k1-x_k))/(z-y);if (fabs(x_k2-x_k1)<EPS){printf("THE ROOT IS x=%g,k=%d\n",x_k2,i);return;}x_k=x_k1;x_k1=x_k2;}printf("AFTER %d repeate,no solved.\n",MAXREPT); } 结果高斯顺序消元法#include<stdio.h>#include<math.h>#define N 4static double aa[N][N+1]={{2,4,0,1,1},{3,8,2,2,3},{1,3,3,0,6},{2,5,2,2,3}}; int gauss(double a[][N+2],double x[]);void putout(double a[][N+2]);void main(){int i,j,det;double a[N+1][N+2],x[N+1];for(i=1;i<=N;i++)for(j=1;j<=N+1;j++)a[i][j]=aa[i-1][j-1];det=gauss(a,x);if(det!=0)for(i=1;i<=N;i++)printf(" x[%d]=%g",i,x[i]);printf("\n");}int gauss(double a[][N+2],double x[]){int i,j,k;double c;putout(a);for(k=1;k<=N-1;k++){ if(fabs(a[k][k])<1e-17){printf("\n pivot element is 0.fail!\n");return 0;}for(i=k+1;i<=N;i++){c=a[i][k]/a[k][k];for(j=k;j<=N+1;j++){a[i][j]=a[i][j]-c*a[k][j];}}putout(a);}if(fabs(a[N][N])<1e-17){printf("\n pivot element is 0.fail!\n");return 0;}for(k=N;k>=1;k--){x[k]=a[k][N+1];for(j=k+1;j<=N;j++){x[k]=x[k]-a[k][j]*x[j];}x[k]=x[k]/a[k][k];}return(1);}void putout(double a[][N+2]){for(int i=1;i<=N;i++){for(int j=1;j<=N+1;j++)printf("%-15g",a[i][j]);printf("\n");}printf("\n");}结果雅克比迭代法#include<stdio.h>#include<math.h>#define N 5#define EPS 0.5e-4static double aa[N][N]={{4,-1,0,-1,0},{-1,4,-1,0,-1},{0,-1,4,-1,0},{-1,0,-1,4,-1},{0,-1,0,-1,4}}; static double bb[N]={2,1,2,1,2};void main(){int i,j,k,NO;double a[N+1][N+1],b[N+1],x[N+1],y[N+1];double d,sum,max;for(i=1;i<=N;i++){for(j=1;j<=N;j++)a[i][j]=aa[i-1][j-1];b[i]=bb[i-1];}printf("\n 请输入最大迭代次数（尽量取大值！）NO:");scanf("%d",&NO);printf("\n");for(i=1;i<=N;i++)x[i]=0;k=0;printf(" k",' ');for(i=1;i<=N;i++)printf("%8cx[%d]",' ',i);printf("\n 0");for(i=1;i<=N;i++)printf("%12.8g",x[i]);printf("\n");do{for(i=1;i<=N;i++){sum=0.0;for(j=1;j<=N;j++)if(j!=i) sum=sum+a[i][j]*x[j];y[i]=(-sum+b[i])/a[i][i];}max=0.0;for(i=0;i<=N;i++){d=fabs(y[i]-x[i]);if(max<d) max=d;x[i]=y[i];}printf("%6d",k+1);for(i=1;i<=N;i++)printf("%12.8g",x[i]);printf("\n");k++;}while((max>=EPS)&&(k<NO));printf("\nk=%d\n",k);if(k>=NO) printf("\nfail!\n");elsefor(i=1;i<=N;i++)printf("x[%d]=%g\t",i,x[i]);}结果拉格朗日插值多项式#include<stdio.h>#include<math.h>#define N 4doublex[N]={0.56160,0.56280,0.56401,0.56521},y[N]={0.82741,0.82659,0.82577,0.82495}; void main(){double x=0.5635;double L(double xx);double lagBasis(int k,double xx);void output();output();printf("\n近似值L(%g)=%g\n",x,L(x));}double lagBasis(int k,double xx){double lb=1;int i;for(i=0;i<N;i++)if(i!=k) lb*=(xx-x[i])/(x[k]-x[i]);return lb;}double L(double xx){double s=0;int i;for(i=0;i<=N;i++)s+=lagBasis(i,xx)*y[i];return s;}void output(){int i;printf("\n各节点信息:\nxi:");for(i=0;i<N;i++)printf("\t%g",x[i]);printf("\nyi:");for(i=0;i<N;i++)printf("\t%g",y[i]);}结果牛顿插值多项式#include <math.h>#include <stdio.h>int a;#define M 4double x[M+1]={0.4,0.55,0.65,0.8,0.9},y[M+1]={0.41075,0.57815,0.69675,0.88811,1.02652}; void main(){double x;printf("输入x=");scanf("%lf",&x);printf("次数：");scanf("%d",&a);double N(double xx,int a);void output();output();printf("\n%d次牛顿插值多项式N(%g)=%g\n",a,x,N(x,a));}double N(double xx,int a){double s=y[0],d=1;int i,j;double df[M+1][M+1];for(i=0;i<=M;i++)df[i][0]=y[i];for(j=1;j<=a;j++)for(i=j;i<=a;i++)df[i][j]=(df[i][j-1]-df[i-1][j-1])/(x[i]-x[i-j]);printf("\nx\tf(x)\t");for(j=1;j<=a;j++) printf("%5d阶",j);for(i=0;i<=a;i++){printf("\n%g\t%g",x[i],y[i]);for(j=1;j<=i;j++)printf("\t%7.5g",df[i][j]);}for(i=1;i<=a;i++){d*=(xx-x[i-1]);s+=df[i][i]*d;}return s;}void output(){int i;printf("\n各节点信息:\nxi:");for(i=0;i<=M;i++)printf("\t%7g",x[i]);printf("\nyi:");for(i=0;i<=M;i++)printf("\t%7g",y[i]);}结果复合梯形公式#include<stdio.h>#include<math.h>#define f(x) 1/(x*x+1)#define Pi 3.1415926void main(){double a=0,b=1;double T,h,x;int n,i;printf("please input n:");scanf("%d",&n);h=(b-a)/n;x=a;T=0;for(i=1;i<n;i++){x+=h;T+=f(x);}T=(f(a)+2*T+f(b))*h/2;printf("T(%d)=%g\n",n,T);printf("The exact value is %g\n",Pi/4);}复合辛普森公式#include<stdio.h>#include<math.h>#define f(x) 1/(1+x*x)#define Pi 3.1415926void main(){double a=0,b=1;double S,h,x;int n,i;printf("please input Even n:");scanf("%d",&n);h=(b-a)/n;x=a; S=0;for(i=1;i<n;i++){x+=h;if(i%2==0) S+=2*f(x);else S+=4*f(x);}S=(f(a)+S+f(b))*h/3;printf("S(%d)=%g\n",n,S);printf("The exact value is %g\n",Pi/4);}龙贝格公式加速#include<stdio.h>#include<math.h>#define f(x) sin(x)/(1+x)#define M 3void main(){double a=0,b=1;double Long(double a,double b);printf("近似值I=%g\n",Long(a,b));}double Long(double a,double b){int n=1,i=1,j=1;double T[M+1][M+1],h,x,sum;h=b-a;T[0][0]=(f(a)+f(b))/2;for(j=1;j<=3;j++){x=a;h/=2;n*=2;sum=0;for(i=1;i<=n;i+=2){x=a+i*h;sum+=f(x);}T[j][0]=T[j-1][0]/2+h*sum;}for(i=1;i<=M;i++)for(j=1;j<=i;j++){T[i][j]=(pow(4,j)*T[i][j-1]-T[i-1][j-1])/(pow(4,j)-1);}printf("k\tT0\tT1\tT2\tT3\n");for(i=0;i<=M;i++){printf("%d",i);for(j=0;j<=i;j++)printf(" %g",T[i][j]);printf("\n");}return T[M][M];}。
西安交通大学概率论上机实验报告
西安交通大学一、试验目的概率论部分1.了解matlab软件的基本命令与操作；2.熟悉matlab用于描述性统计的基本菜单操作及命令；3.会用matlab求密度函数值、分布函数值、随机变量分布的上下侧分位数。
数理统计部分1.熟悉matlab进行参数估计、假设检验的基本命令与操作.2.掌握用matlab生成点估计量值的模拟方法3.会用matlab进行总体数学期望和方差的区间估计。
4.会用matlab进行单个、两个正态总体均值的假设检验。
5.会用matlab进行单个、两个正态总体方差的假设检验。
二、试验问题实验五、随机变量综合试验实验内容1. 产生(6),(10),F(6,10)和t(6)四种随机数，并画出相应的频率直方图；2. 在同一张图中画出了N(0,1)和t(6)随机数频率直方图，比较它们的异同；3. 写出计算上述四种分布的分布函数值和相应上侧分位点命令.实验七、对统计中参数估计进行计算机模拟验证实验内容：1.产生服从给定分布的随机数，模拟密度函数或概率分布；2.对分布包含的参数进行点估计，比较估计值与真值的误差；3. 对分布包含的参数进行区间估计，行区间估计，可信度。
三、实验源程序及结果实验5源程序：% 清空内存，清空输出屏幕clc;clear;% 首先是指数分布n = normpdf(-2::14,6);% 绘制频率直方图plot(-2::14,n,'color','r','linewidth',2);ylabel('概率密度');title('正态分布概率密度');% t分布h1 = figure;t = tpdf(-3::3,6);plot(-3::3,t,'color','g','linewidth',2);ylabel('对应频率');title('t分布频率密度');% F分布h2 = figure;f = fpdf(0::10,6,10);plot(0::10,f,'color','k','linewidth',2); ylabel('对应频率');title('F分布频率直方图');% 卡方分布h3 = figure;ka = chi2pdf(0::15,6);plot(0::15,ka,'color','y','linewidth',2); ylabel('对应频率');title('卡方分布频率直方图');% 再来绘图h4 = subplot(2,1,1);y1=normpdf(-10::10,0,1);plot(-10::10,y1,'color','b','linewidth',2); title('N(0,1)');h5 = subplot(2,1,2);t1 = tpdf(-10::10,6);plot(-10::10,t1,'color','r','linewidth',2);%上侧分位数norminv,0,1)tinv,6)chi2inv,6)finv,6,10)运行结果：正态分布T分布F分布N(0,1)和t(6)随机数频率直方图四种分布的分布函数值和相应上侧分位点实验7源程序：% 以正太分布为例% 清空内存，清空输出屏幕clc;clear;y=normrnd(10,1,10000,1);ymin=min(y);ymax=max(y);x=linspace(ymin,ymax,80);yy=hist(y,x);yy=yy/10000;bar(x,yy);grid;xlabel('(a)¸概率密度分布直方图 ');phat=mle(y,'distribution','norm','alpha',%对分布函数参数进行区间估计，并估计区间的可信度 [mu,sigma,m_ci,s_si]=normfit(y,运行结果：正态分布概率密度分布直方图得到估计参数m =σ =由上可知估计的m = ，而实际是 10。
西安交大概率论上机实验报告西安交通大学概率论实验报告
西安交⼤概率论上机实验报告西安交通⼤学概率论实验报告概率论与数理统计上机实验报告⼀、实验内容使⽤MATLAB 软件进⾏验证性实验，掌握⽤MATLAB 实现概率统计中的常见计算。
本次实验包括了对⼆维随机变量，各种分布函数及其图像以及频率直⽅图的考察。
1、列出常见分布的概率密度及分布函数的命令，并操作。
2、掷硬币150次，其中正⾯出现的概率为0.5，这150次中正⾯出现的次数记为X ，(1) 试计算45=X 的概率和45≤X 的概率；(2) 绘制分布函数图形和概率分布律图形。
3、⽤Matlab 软件⽣成服从⼆项分布的随机数，并验证泊松定理。
4、设22221),(y x e y x f +-=π是⼀个⼆维随机变量的联合概率密度函数，画出这⼀函数的联合概率密度图像。
5、来⾃某个总体的样本观察值如下，计算样本的样本均值、样本⽅差、画出频率直⽅图。
A=[16 25 19 20 25 33 24 23 20 24 25 17 15 21 22 26 15 23 2220 14 16 11 14 28 18 13 27 31 25 24 16 19 23 26 17 14 30 21 18 16 18 19 20 22 19 22 18 26 26 13 21 13 11 19 23 18 24 28 13 11 25 15 17 18 22 16 13 12 13 11 09 15 18 21 15 12 17 13 14 12 16 10 08 23 18 11 16 28 13 21 22 12 08 15 21 18 16 16 19 28 19 12 14 19 28 28 28 13 21 28 19 11 15 18 24 18 16 28 19 15 13 22 14 16 24 20 28 18 18 28 14 13 28 29 24 28 14 18 18 18 08 21 16 24 32 16 28 19 15 18 18 10 12 16 26 18 19 33 08 11 18 27 23 11 22 22 13 28 14 22 18 26 18 16 32 27 25 24 17 17 28 33 16 20 28 32 19 23 18 28 15 24 28 29 16 17 19 18]6. 利⽤Matlab 软件模拟⾼尔顿板钉试验。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
计算方法上机报告姓名：学号：班级：能动上课班级：题目及求解：一、对以下和式计算：∑∞⎪⎭⎫⎝⎛+-+-+-+=0681581482184161n n n n S n ，要求：① 若只需保留11个有效数字，该如何进行计算； ② 若要保留30个有效数字，则又将如何进行计算；1 算法思想（1）根据精度要求估计所加的项数，可以使用后验误差估计，通项为： 142111416818485861681n n na n n n n n ε⎛⎫=---<< ⎪+++++⎝⎭； （2）为了保证计算结果的准确性，写程序时，从后向前计算； （3）使用Matlab 时，可以使用以下函数控制位数： digits(位数)或vpa(变量，精度为数)2 算法结构;0=s⎪⎭⎫⎝⎛+-+-+-+=681581482184161n n n n t n ； for 0,1,2,,n i =⋅⋅⋅ if 10m t -≤end;for ,1,2,,0n i i i =--⋅⋅⋅;s s t =+3 Matlab 源程序clear; %清除工作空间变量 clc; %清除命令窗口命令 m=input('请输入有效数字的位数m='); %输入有效数字的位数 s=0;for n=0:50t=(1/16^n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6));if t<=10^(-m) %判断通项与精度的关系break;endend;fprintf('需要将n值加到n=%d\n',n-1); %需要将n值加到的数值for i=n-1:-1:0t=(1/16^i)*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6));s=s+t; %求和运算ends=vpa(s,m) %控制s的精度4 结果与分析若保留11位有效数字，则n=7，此时求解得：s =3.1415926536；若保留30位有效数字时，则n=22, 此时求解得：s =3.8。
通过上面的实验结果可以看出，通过从后往前计算，这种算法很好的保证了计算结果要求保留的准确数字位数的要求。
二、某通信公司在一次施工中，需要在水面宽度为20米的河沟底部沿直线走向铺设一条沟底光缆。
在铺设光缆之前需要对沟底的地形进行初步探测，从而估计所需光缆的长度，为工程预算提供依据。
已探测到一组等分点位置的深度数据(单位:米)如下表所示：分点0 1 2 3 4 5 6深度9.01 8.96 7.96 7.97 8.02 9.05 10.13分点7 8 9 10 11 12 13深度11.18 12.26 13.28 13.32 12.61 11.29 10.22分点14 15 16 17 18 19 20深度9.15 7.90 7.95 8.86 9.81 10.80 10.93①请用合适的曲线拟合所测数据点；②预测所需光缆长度的近似值，作出铺设河底光缆的曲线图；1 算法思想如果使用多项式差值，则由于龙格现象，误差较大，因此，用相对较少的插值数据点作插值，可以避免大的误差，但是如果又希望将所得数据点都用上，且所用数据点越多越好，可以采用分段插值方式，即用分段多项式代替单个多项式作插值。
分段多项式是由一些在相互连接的区间上的不同多项式连接而成的一条连续曲线，其中三次样条插值方法是一种具有较好“光滑性”的分段插值方法。
在本题中，假设所铺设的光缆足够柔软，在铺设过程中光缆触地走势光滑，紧贴地面，并且忽略水流对光缆的冲击。
海底光缆线的长度预测模型如图2-1所示，光缆从Ah。
点铺至B点，在某点处的深度为i图2-1 海底光缆线的长度预测模型计算光缆长度时，用如下公式：20()L f x ds =⎰200(f x =⎰ 191(k kk f x +==∑⎰=2 算法结构1） For n i ,,2,1,0⋅⋅⋅=1.1 i i M y ⇒ 2） For 2,1=k2.1 For k n n i ,,1, -=2.1.1 i k i i i i M x x M M ⇒----)/()(13）101h x x ⇒-4）For 1-,,2,1n i =4.1 11++⇒-i i i h x x4.2 b a c c h h h i i i i i i ⇒⇒-⇒+++2;1;)/(11 4.3 i i d M ⇒+165）0000;;c M d M d n n ⇒⇒⇒λn n n b a b ⇒⇒⇒2;;20μ6）1111,γμ⇒⇒d b7）获取M 的矩阵元素个数，存入m 8）For m k ,,3,2 =8.1 k k k l a ⇒-1/μ8.2 k k k k c l b μ⇒⋅-1- 8.3 k k k k l d γγ⇒⋅-1- 9）m m m M ⇒μγ/10）For 1,,2,1 --=m m k10.1 k k k k k M M c ⇒⋅-+μγ/)(1 11）获取x 的元素个数存入s 12）k ⇒113） For 1,,2,1-=s i13.1 if i x x ≤~then k i ⇒；break else k i ⇒+114）xx x x x x h x x k k k k ˆ~;~;11⇒-⇒-⇒--- y h x h M y x h M y x M x M k k k k k k ~/]ˆ)6()6(6ˆ6[2211331⇒-+-++---3 Matlab 源程序clear; %清除工作空间变量 clc; %清除命令窗口命令x=0:1:20; %产生从0到20含21个等分点的数组 X=0:0.2:20;y=[9.01,8.96,7.96,7.97,8.02,9.05,10.13,11.18,12.26,13.28,13.32,12.61,11.29,10.22,9.15,7.90,7.95,8.86,9.81,10.80,10.93]; %等分点位置的深度数据 n=length(x); %等分点的数目 N=length(X);%% 求三次样条插值函数s(x) M=y;for k=2:3; %计算二阶差商并存放在M 中 for i=n:-1:k;M(i)=(M(i)-M(i-1))/(x(i)-x(i-k+1)); end endh(1)=x(2)-x(1); %计算三对角阵系数a,b,c及右端向量d for i=2:n-1;h(i)=x(i+1)-x(i);c(i)=h(i)/(h(i)+h(i-1));a(i)=1-c(i);b(i)=2;d(i)=6*M(i+1);endM(1)=0; %选择自然边界条件M(n)=0;b(1)=2;b(n)=2;c(1)=0;a(n)=0;d(1)=0;d(n)=0;u(1)=b(1); %对三对角阵进行LU分解y1(1)=d(1);for k=2:n;l(k)=a(k)/u(k-1);u(k)=b(k)-l(k)*c(k-1);y1(k)=d(k)-l(k)*y1(k-1);endM(n)=y1(n)/u(n); %追赶法求解样条参数M(i)for k=n-1:-1:1;M(k)=(y1(k)-c(k)*M(k+1))/u(k);ends=zeros(1,N);for m=1:N;k=1;for i=2:n-1if X(m)<=x(i);k=i-1;break;elsek=i;endendH=x(k+1)-x(k); %在各区间用三次样条插值函数计算X点处的值x1=x(k+1)-X(m);x2=X(m)-x(k);s(m)=(M(k)*(x1^3)/6+M(k+1)*(x2^3)/6+(y(k)-(M(k)*(H^2)/6))*x1+(y(k+1)-(M(k+1)*(H^2)/6 ))*x2)/H;end%% 计算所需光缆长度L=0; %计算所需光缆长度for i=2:NL=L+sqrt((X(i)-X(i-1))^2+(s(i)-s(i-1))^2);enddisp('所需光缆长度为L=');disp(L);figureplot(x,y,'*',X,s,'-') %绘制铺设河底光缆的曲线图xlabel('位置','fontsize',16); %标注坐标轴含义ylabel('深度/m','fontsize',16);title('铺设河底光缆的曲线图','fontsize',16);grid;4 结果与分析铺设海底光缆的曲线图如图2-2所示：图2-2 铺设河底光缆曲线图仿真结果表明，运用分段三次样条插值所得的拟合曲线能较准确地反映铺设光缆的走势图，计算出所需光缆的长度为 L=26.4844m 。
三、 假定某天的气温变化记录如下表所示，试用数据拟合的方法找出这一天的气温变化的规律；试计算这一天的平均气温，并试估计误差。
1 算法思想在本题中，数据点的数目较多。
当数据点的数目很多时，用“多项式插值”方法做数据近似要用较高次的多项式，这不仅给计算带来困难，更主要的缺点是误差很大。
用“插值样条函数”做数据近似，虽然有很好的数值性质，且计算量也不大，但存放参数Mi 的量很大，且没有一个统一的数学公式来表示，也带来了一些不便。
另一方面，在有的实际问题中，用插值方法并不合适。
当数据点的数目很大时，要求)(x p 通过所有数据点，可能会失去原数据所表示的规律。
如果数据点是由测量而来的，必然带有误差，插值法要求准确通过这些不准确的数据点是不合适的。
在这种情况下，不用插值标准而用其他近似标准更加合理。
通常情况下，是选取i α使2E 最小，这就是最小二乘近似问题。
在本题中，采用“最小二乘法”找出这一天的气温变化的规律，使用二次函数、三次函数、四次函数以及指数型函数2)(c t b ae C --=，计算相应的系数，估算误差，并作图比较各种函数之间的区别。
2 算法结构本算法用正交化方法求数据的最小二乘近似。
假定数据以用来生成了G ，并将y 作为其最后一列（第1+n 列）存放。
结果在a 中，ρ是误差22E 。
I 、使用二次函数、三次函数、四次函数拟合时 1. 将“时刻值”存入i x ，数据点的个数存入m2. 输入拟合多项式函数)(x p 的最高项次数1-=n h ，则拟合多项式函数为)()()()(2211x g x g x g x p n n ααα+⋅⋅⋅++= , 根据给定数据点确定G For1,,2,1,0-⋅⋅⋅=n jFor m i ,,2,1⋅⋅⋅=2.1 1,+⇒j i ji g x 2.2 1,+⇒n i i g y 3. For n k ,,2,1⋅⋅⋅=3.1 [形成矩阵k Q ]3.1.1 σ⇒-∑=mk i ikkk g g 2/12))(sgn( 3.1.2 k kk g ωσ⇒-3.1.3 For m k k j ,,2,1⋅⋅⋅++=3.1.3.1 j jk g ω⇒ 3.1.4 βσω⇒k 3.2 [变换1-k G 到k G ] 3.2.1 kk g ⇒σFor 1,,,2,1+⋅⋅⋅++=n n k k j 3.2.2 t g mk i ij i ⇒∑=βω/)(3.2.3 For m k k i ,,1,⋅⋅⋅+=3.2.3.1 ij i ij g t g ⇒+ω4. [解三角方程1h Ra =] 4.1nnn n n a g g ⇒+/1.4.2 For 1,,2,1⋅⋅⋅--=n n i 4.2.1iii ni j j ijn i a g x gg ⇒-∑+=+/][11.5. [计算误差22E ]ρ⇒∑+=+mn i n i g121,II 、使用指数函数拟合时现将指数函数进行变形： 将y C =,x t=代入2)(c t b aeC --=得:2)(c x b aey --=对上式左右取对数得: 222ln ln bx bcx bc a y -+-=令b bc bc a y z -==-==21202ln ln ααα，，，则可得多项式: 2210x x z ααα++=（3）Matlab 源程序clear; %清除工作空间变量 clc; %清除命令窗口命令 x=0:24; %将时刻值存入数组 y=[15,14,14,14,14,15,16,18,20,20,23,25,28,31,34,31,29,27,25,24,22,20,18,17,16]; [~,m]=size(x); %将数据点的个数存入m T=sum(y(1:m))/m;fprintf('一天的平均气温为T=%f\n',T); %求一天的平均气温 %% 二次、三次、四次函数的最小二乘近似h=input('请输入拟合多项式的最高项次数='); %根据给定数据点生成矩阵G n=h+1; G=[];for j=0:(n-1)g=x.^j; %g(x)按列排列G=vertcat(G,g); %g垂直连接GendG=G'; %转置得到矩阵Gfor i=1:m %将数据y作为G的最后一列（n+1列）G(i,n+1)=y(i);endG;for k=1:n %形成矩阵Q(k)if G(k,k)>0;sgn=1;elseif G(k,k)==0;sgn=0;else sgn=-1;endsgm=-sgn*sqrt(sum(G(k:m,k).^2));w=zeros(1,n);w(k)=G(k,k)-sgm;for j=k+1:mw(j)=G(j,k);endbt=sgm*w(k);G(k,k)=sgm; %变换Gk-1到Gkfor j=k+1:n+1t=sum(w(k:m)*G(k:m,j))/bt;for i=k:m;G(i,j)=G(i,j)+t*w(i);endendendA (n)=G(n,n+1)/G(n,n); %解三角方程求系数Afor i=n-1:-1:1A (i)=(G(i,n+1)-sum(G(i,i+1:n).*A (i+1:n)))/G(i,i);ende=sum(G(n+1:m,n+1).^2); %计算误差efprintf('%d次函数的系数是：',h); %输出系数a及误差edisp(A);fprintf('使用%d次函数拟合的误差是%f：',h,e);t=0:0.05:24;A=fliplr(A); %将系数数组左右翻转Y=poly2sym(A); %将系数数组转化为多项式subs(Y,'x',t);Y=double(ans);figure(1)plot(x,y,'k*',t,Y,'r-'); %绘制拟合多项式函数图形xlabel('时刻'); %标注坐标轴含义ylabel('平均气温');title([num2str(n-1),'次函数的最小二乘曲线']);grid;%% 指数函数的最小二乘近似yy=log(y);n=3;G=[];GG=[];for j=0:(n-1)g=x.^j; %g(x)按列排列G=vertcat(G,g); %g垂直连接Ggg=t.^j; %g(x)按列排列GG=vertcat(GG,gg); %g垂直连接GendG=G'; %转置得到矩阵Gfor i=1:m %将数据y作为G的最后一列（n+1列）G(i,n+1)=yy(i);endG;for k=1:n %形成矩阵Q(k)if G(k,k)>0;sgn=1;elseif G(k,k)==0;sgn=0;else sgn=-1;endsgm=-sgn*sqrt(sum(G(k:m,k).^2));w=zeros(1,n);w(k)=G(k,k)-sgm;for j=k+1:mw(j)=G(j,k);endbt=sgm*w(k);G(k,k)=sgm; %变换Gk-1到Gkfor j=k+1:n+1t=sum(w(k:m)*G(k:m,j))/bt;for i=k:m;G(i,j)=G(i,j)+t*w(i);endendendA(n)=G(n,n+1)/G(n,n); %解三角方程求系数A for i=n-1:-1:1A (i)=(G(i,n+1)-sum(G(i,i+1:n).*A (i+1:n)))/G(i,i);endb=-A(3);c=A(2)/(2*b);a=exp(A(1)+b*(c^2));G(n+1:m,n+1)=exp(sum(G(n+1:m,n+1).^2));e=sum(G(n+1:m,n+1).^2); %计算误差efprintf('\n指数函数的系数是：a=%f,b=%f,c=%f',a,b,c); %输出系数及误差e fprintf('\n使用指数函数拟合的误差是：%f',e);t=0:0.05:24;YY=a.*exp(-b.*(t-c).^2);figure(2)plot(x,y,'k*',t,YY,'r-'); %绘制拟合指数函数图形xlabel('时刻'); %标注坐标轴含义ylabel('平均气温');title(['指数函数的最小二乘曲线']);grid;4 结果与分析a 二次函数：一天的平均气温为：21.2000二次函数的系数: 8.3063 2.6064 -0.0938使用二次函数拟合的误差是：280.339547二次函数的最小二乘曲线如下图3-1所示：图3-1 二次函数的最小二乘曲线b 三次函数：一天的平均气温为：21.2000三次函数的系数: 13.3880 -0.2273 0.2075 -0.0084 使用三次函数拟合的误差是：131.061822三次函数的最小二乘曲线如图3-2所示：图3-2 三次函数的最小二乘曲线c四次函数：一天的平均气温为：21.2000四次函数的系数: 16.7939 -3.7050 0.8909 -0.0532 0.0009 使用四次函数拟合的误差是：59.04118四次函数的最小二乘曲线如图3-3所示：图3-3 四次函数的最小二乘曲线d 指数函数：一天的平均气温为：21.2000指数函数的系数是： a=26.160286，b=0.004442，c=14.081900 使用指数函数拟合的误差是： 57.034644 指数函数的最小二乘曲线如下图所示：图3-4 指数函数的最小二乘曲线通过上述几种拟合可以发现，多项式的次数越高，计算拟合的效果越好，误差越小，说明结果越准确；同时，指数多项式拟合的次数虽然不高，但误差最小，说明结果最准确。