计算方法上机作业集合Word格式文档下载.docx

1、上机作业181页第四题线性方程组为（1）顺序消元法A=1.1348,3.8326,1.1651,3.4017;0.5301,1.7875,2.5330,1.5435;3.4129,4.9317,8.7643,1.3142;1.2371,4.9998,10.6721,0.0147;b=9.5342;6.3941;18.4231;16.9237;上机代码函数部分如下function b = gaus（ A,b ）%用b返回方程组的解B=A,b;n=length（b）;RA=rank（A）;RB=rank（B）;dif=RB-RA;if dif disp（此方程组无解）; returnendif R

2、A=RB if RA=n format long;此方程组有唯一解 for p=1:n-1 for k=p+1:n m=B（k,p）/B（p,p）; B（k,p:n+1）=B（k,p:n+1）-m*B（p,p:n+1）; end end %顺序消元形成上三角矩阵 b=B（1:n,n+1）; A=B（1:n,1:n）; b（n）=b（n）/A（n,n）; for q=n-1:-1:1 b（q）=（b（q）-sum（A（q,q+1:n）*b（q+1:n）/A（q,q）; end %回代求解 else此方程组有无数组解上机运行结果为 A=1.1348,3.8326,1.1651,3.4017; X=

3、gaus（A,b）此方程组有唯一解X = 1.000000000000000（2）列主元消元法函数部分代码如下function b = lieZhu（A,b ）%用b返回方程组的解format long;该方程组无解if dif=0该方程组有唯一解 c=zeros（1,n）; for i=1: max=abs（A（i,i）; m=i; for j=i+1: if max0.0001）%判断当前的解是否到达精度要求 b=x1; x1=B*b+f; num=num+1;end;fprintf（求得的解为：ndisp（x1）;迭代次数：%d次n,num）;上机运行，结果如下 A=4,-1,0,-1,

4、0,0;-1,4,-1,0,-1,0;0,-1,4,-1,0,-1;-1,0,-1,4,-1,0;0,-1,0,-1,4,-1;0,0,-1,0,-1,4; b=0;5;-2;6; Jacobi（A,b） 0.999981765074381 1.99995018125674 0.999975090628368 1.9999635301487628次满足要求的解如输出结果所示，总共迭代了28次（2）高斯-赛德尔迭代法上机程序如下所示function =Gauss_Seidel（ A,b ）%高斯赛德尔迭代法D=diag（diag（A）;L=D-tril（A）;U=D-triu（A）;B=inv（

5、D-L）*U;f=inv（D-L）*b;x0=B*b+f;while（norm（x0-b）0.0001） b=x0; x0=B*b+f;结果为ndisp（x0）;迭代次数为： Gauss_Seidel（A,b）结果为 0.999960143810658 1.99995676152139 0.999963508299833 1.99996662162874 0.999972527179715 1.9999840088698915次满足精度要求的解如上述程序打印输出所示，迭代了15次（3）SOR迭代法w依次取1.334,1.95,0.95上机代码如下function = SOR（A,b,w ）%S

6、OR迭代法B=inv（D-w*L）*（1-w）*D+w*U）;f=w*inv（D-w*L）*b;迭代次数为%dn 上机运行 SOR（A,b,1.334） 1.00001878481009 1.99998698322858 1.00001815013068 2.00000041318053 0.999991858543476 2.0000068413569迭代次数为13 SOR（A,b,1.95） 0.999984952088107 2.00000960832604 0.999959115182729 2.0000168426006 1.00000443526697 1.999978851134

7、46迭代次数为241 SOR（A,b,0.95） 0.999961518309351 1.99995706825231 0.999963054838453 1.99996580572033 0.999971141727589 1.9999827300678迭代次数为17由以上输出得到w取值不同的情况下，得到的满足精度要求的结果，迭代次数分别如输出所示第五次上机作业8.从函数表x0.00.10.20.30.4010.5f（x）0.398940.396950.391420.381380.368120.35206出发，用以下方法计算f（0.15）,f（0.31）及f（0.47）的近似值（1）分段线性

8、插值（2）分段二次插值（3）全区间上拉格朗日插值1线性插值编写函数如下function R = div_line（ x0,y0,x ）%线性插值p=length（x0）;n=length（y0）;m=length（x）;if（p=n）%x的个数与y的个数不等 error（数据输入有误，请重新输入 return; fprintf（线性插值n for t=1:m z=x（t）; if（zx0（p）x%d不在所给自变量范围内，无法进行插值,t）; continue;p-1x0（i+1） break; end; R（t）=y0（i）*（x（t）-x0（i+1）/（x0（i）-x0（i+1）+y0（i+

9、1）*（x（t）-x0（i）/（x0（i+1）-x0（i）;上机运行如下 x0=0.0 0.1 0.195 0.3 0.401 0.5; y0=0.39894 0.39695 0.39142 0.38138 0.36812 0.35206; x=0.15 0.31 0.47; div_line（x0,y0,x）线性插值ans = 0.39404 0.38007 0.35693即结果为f（0.15）0.39404,f（0.31）0.38007,f（0.47）0.356932分段二次插值编写的函数如下function R = div2line（x0,y0,x）%分段二次插值m=length（y0）

10、;n=length（x）;if（p=m）输入错误，请重新输入数据for t=1: if（x（t）x%d不在所给区间上 tag=2; m=abs（x（t）-x0（1）+abs（x（t）-x0（2）+abs（x（t）-x0（3）; for i=3: sum=abs（x（t）-x0（i-1）+abs（x（t）-x0（i）+abs（x（t）-x0（i+1）; if（summ） m=sum; tag=i;tag=%dn,tag）;R（t）=y0（tag-1）*（x（t）-x0（tag）*（x（t）-x0（tag+1）/（x0（tag-1）-x0（tag）*（x0（tag-1）-x0（tag+1）+y0（

11、tag）*（x（t）-x0（tag-1）*（x（t）-x0（tag+1）/（x0（tag）-x0（tag-1）*（x0（tag）-x0（tag+1）+y0（tag+1）*（x（t）-x0（tag-1）*（x（t）-x0（tag）/（x0（tag+1）-x0（tag-1）*（x0（tag+1）-x0（tag）;End上机运行，执行结果为div2line（x0,y0,x） 0.39448 0.38022 0.35725即分段二次插值方法下，f（0.15）0.39448,f（0.31）0.38022,f（0.47）0.357253上机编写的程序如下function R = lagrange（x0,y

12、0,x）%全区间上拉格朗日插值p=length（y0）;n=length（x0）;%计算函数表和x的长度if p = n error（%假设函数表的x与y长度不一致则输入有误else fprintf（拉格朗日插值nm %利用循环计算每个x的插值 s=0.0; for k=1:n p=1; if i=k p=p*（z-x0（i）/（x0（k）-x0（i）; s=s+y0（k）*p; %根据拉格朗日插值公式求解yR（t）=s;%输出插值结果 lagrange（x0,y0,x）拉格朗日插值 0.39447 0.38022 0.357220.39447,f（0.31）0.357229.解：上机程序如下

13、，为方便起见，将所有操作分在四个函数中进行入口函数function =spline（ X,Y,xx,y1_0,y1_18 ）%输出自变量所对应的函数值M=getM（X,Y,y1_0,y1_18）;%先得到Ms=xx;k=length（xx）;for a=1:k s（xx（a）=getExp（X,Y,M,xx（a）;%计算自变量所在小区间对应曲线的表达式，并根据表达式计算对应的函数值s（%d）=%fn,xx（a）,s（xx（a）; %输出打印函数值获取Mfunction M = getM（X,Y,y1_0,y1_1）%得到Mn=length（X）;a=0*X;b=a;c=a;h=a;f=a;b=

14、b+2;h（2:n）=X（2:n）-X（1:n-1）; % h（1）不用 a（2:n-1）=h（2:n-1）./（h（2:n-1）+h（3:n）;c（2:n-1）=1-a（2:a（n）=1;c（1）=1;Yf（2:n）=Y（2:n）-Y（1:f（2:n-1）=6*（Yf（3:n）./h（3:n）-Yf（2:n-1）./h（2:n-1）./（h（2:f（1）=6*（Yf（2）/h（2）-y1_0）/h（2）;f（n）=6*（y1_1-Yf（n）/h（n）/h（n）;M=CalM（n,a,b,c,f）;%计算M计算Mfunction f = CalM（n,a,b,c,f）% 追赶法求解Meps=0

15、.1e-15; %防止参数过小，是的计算误差过大if abs（b（1）eps除数为0，停止计算 c（1）=c（1）/b（1）;%追赶法：根据递推算式计算for i=2: b（i）=b（i）-a（i）*c（i-1）; if abs（b（i）X（n5） n1=n5; n2=n5;%计算yy=M（n1）*（X（n2）-x）3/（6*h（n2）+ M（n2）*（x-X（n1）3/（6*h（n2）;y=y+（Y（n1）-M（n1）*h（n2）*h（n2）/6）*（X（n2）-x）/h（n2）;y=y+（Y（n2）-M（n2）*h（n2）*h（n2）/6）*（x-X（n1）/h（n2）;%结束上机试运行，

16、过程如下 X=0.52 3.1 8.0 17.95 28.65 39.62 50.65 78 104.6 156.6 208.6 260.7 312.5 364.4 416.3 468 494 507 520; Y=5.28794 9.4 13.84 20.2 24.9 28.44 31.1 35 36.5 36.6 34.6 31.0 26.34 20.9 14.8 7.8 3.7 1.5 0.2; xx=2 4 6 12 16 30 60 110 180 280 400 515; y1_0=1.86548; y1_18=-0.046115;spline（X,Y,xx,y1_0,y1_18）s（2）=7.825123s（4）=10.481311s（6）=12.363477s（12）=16.575574s（16）=19.091594s（30）=25.386402s（60）=32.804283s（110）=36.647878s（180）=35.917148s（280）=29.368462s（400）=16.799784s（515）=0.542713根据上述程序运行结果，可得所有自变量对应的函数值，如上输出结果所示第六次上机作业10.已知一组实验数据i23456789