2022年2022年计算方法上机作业集合 .pdf

《2022年2022年计算方法上机作业集合 .pdf》由会员分享，可在线阅读，更多相关《2022年2022年计算方法上机作业集合 .pdf（29页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、第一次 &第二次上机作业上机作业：1.在 Matlab 上执行： 5.1-5-0.1 和 1.5-1-0.5 给出执行结果，并简要分析一下产生现象的原因。解：执行结果如下：在 Matlab 中，小数值很难用二进制进行描述。由于计算精度的影响，相近两数相减会出现误差。2.（课本 181 页第一题）解： （1）n=0 时，积分得 ?0=ln6-ln5，编写如下图代码名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 1 页，共 29 页 - - - - - - - - - 从 以 上 代 码

2、 显 示 的 结 果 可 以 看 出 ， ?20的 近 似 值 为0.012712966517465 （2）?=?5+ ?10dx,可得?610dx?5+ ?10dx?510dx,得16( ? +1)?15(? +1),则有1126?201105, 取?20=1105,以此逆序估算 ?0。代码段及结果如下图所示名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 2 页，共 29 页 - - - - - - - - - 结 果 是 从 ?19逆 序 输 出 至 ?0， 所 以 得 到 ?0

3、的 近 似 值 为0.088392216030227。（3）从 ?20估计的过程更为可靠。首先根据积分得表达式是可知，被积函数随着n 的增大，其所围面积应当是逐步减小的，即积分值应是随着n 的递增二单调减小的， （1）中输出的值不满足这一条件，（2）满足。设?表示 ?的近似值，?-?=(- 5)?(?0- ?0)(根据递推公式可以导出此式),可以看出，随着 n 的增大，误差也在增大，所以顺序估计时，算法不稳定性逐渐增大，逆序估计情况则刚好相反，误差不断减小，算法逐渐趋于稳定。2.（课本 181 页第二题）名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - -

4、- - - - - - 名师精心整理 - - - - - - - 第 3 页，共 29 页 - - - - - - - - - （1）上机代码如图所示求得近似根为0.09058 （2）上机代码如图所示得近似根为0.09064；名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 4 页，共 29 页 - - - - - - - - - （3）牛顿法上机代码如下计算所得近似解为0.09091 第三次上机作业上机作业 181 页第四题线性方程组为1.13483.83260.53011.7875

5、1.16513.40172.53301.54353.41294.93171.23714.99988.76431.314210.67210.0147? 1? 2? 3? 4=9.53426.394118.423116.9237（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(

6、A,b )% 用b返回方程组的解B=A,b; n=length(b); RA=rank(A); 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 5 页，共 29 页 - - - - - - - - - RB=rank(B); dif=RB-RA; ifdif0 disp(此方程组无解 ); return end if RA=RB if RA=n format long; disp(此方程组有唯一解 ); for p=1:n-1 for k=p+1:n m=B(k,p)/B(p,p);

7、 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 disp(此方程组有无数组解 ); end end 上机运行结果为 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.

8、2371,4.9998,10.6721,0.0147; 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 6 页，共 29 页 - - - - - - - - - b=9.5342;6.3941;18.4231;16.9237; X=gaus(A,b) 此方程组有唯一解X = 1.000000000000000 1.000000000000000 1.000000000000000 1.000000000000000 （2）列主元消元法A=1.1348,3.8326,1.1651,3.

9、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 = lieZhu(A,b )% 用 b 返回方程组的解B=A,b; RA=rank(A); RB=rank(B); n=length(b); dif=RB-RA; format long; 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整

10、理 - - - - - - - 第 7 页，共 29 页 - - - - - - - - - ifdif0 disp(该方程组无解 ); return end ifdif=0 if RA=n disp(该方程组有唯一解 ); c=zeros(1,n); fori=1:n-1 max=abs(A(i,i); m=i; for j=i+1:n if max 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.53

11、42;6.3941;18.4231;16.9237; X=lieZhu(A,b) 该方程组有唯一解X = 1.000000000000000 1.000000000000002 0.999999999999999 0.999999999999999 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 9 页，共 29 页 - - - - - - - - - 根据两种方法运算结果， 顺序消元法得到结果比列主元消元法得到的结果精度更高。（注：matlab 使用的是 2015b 版本，不知道是

12、保留小数位数太少，还是程序原因，顺序消元输出结果总是等于准确解，请老师指正）第四次上机作业7. 分析用下列迭代法解线性方程组 4- 1- 1 4 0- 1- 1 0 0 0- 1 0 0- 1- 1 0 4- 1- 1 4 0- 1- 1 00- 10 0 0- 1- 1 0 4- 1- 1 4?1?2?3?4?5?6= 0 5- 2 5- 2 6的收敛性，并求出使 ?(? +1)- ?(? )2 0.0001 的近似解及相应的迭代次数。（1）雅可比迭代法解：上机编写的函数如下function = Jacobi(X,b)% 雅可比迭代法D=diag(diag(X);% 得到对角线元素组成的矩阵

13、B=inv(D)*(D-X);f=inv(D)*b;b(:,:)=0;x1=B*b+f;num=1;while(norm(x1-b)0.0001)% 判断当前的解是否达到精度要求 b=x1;x1=B*b+f;num=num+1;end ; fprintf( 求得的解为： n);disp(x1);名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 10 页，共 29 页 - - - - - - - - - fprintf( 迭代次数： %d次n,num);end上机运行，结果如下 A=4,

14、-1,0,-1,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;5;-2;6; Jacobi(A,b) 求得的解为： 0.999981765074381 1.99995018125674 0.999975090628368 1.99995018125674 0.999975090628368 1.99996353014876 迭代次数： 28 次满足要求的解如输出结果所示，总共迭代了28 次（2）高斯-赛德尔迭代法上机程序如下所示function =Gauss_Sei

15、del( A,b )% 高斯赛德尔迭代法D=diag(diag(A);L=D-tril(A);U=D-triu(A);名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 11 页，共 29 页 - - - - - - - - - B=inv(D-L)*U;f=inv(D-L)*b;b(:,:)=0;x0=B*b+f;num=1;while(norm(x0-b)0.0001)num=num+1; b=x0;x0=B*b+f;end ;fprintf( 结果为 n);disp(x0);fpr

16、intf( 迭代次数为： %d次n,num);end A=4,-1,0,-1,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;5;-2;6; Gauss_Seidel(A,b) 结果为 0.999960143810658 1.99995676152139 0.999963508299833 1.99996662162874 0.999972527179715 1.99998400886989 迭代次数为： 15次名师资料总结 - - -精品资料欢迎下载 - - - -

17、 - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 12 页，共 29 页 - - - - - - - - - 满足精度要求的解如上述程序打印输出所示，迭代了15 次（3）SOR迭代法（ w依次取 1.334,1.95,0.95）上机代码如下function = SOR(A,b,w )%SOR 迭代法 D=diag(diag(A);L=D-tril(A);U=D-triu(A);B=inv(D-w*L)*(1-w)*D+w*U);f=w*inv(D-w*L)*b;b(:,:)=0;x0=B*b+f;num=1;while(norm(x0-b)

18、0.0001)num=num+1; b=x0;x0=B*b+f;end ;fprintf( 结果为 n);disp(x0);fprintf( 迭代次数为 %dn ,num);end上机运行 A=4,-1,0,-1,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;5;-2;6; SOR(A,b,1.334) 结果为名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第

19、 13 页，共 29 页 - - - - - - - - - 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.99997885113446 迭代次数为 241 SOR(A,b,0.95) 结果为 0.9999615

20、18309351 1.99995706825231 0.999963054838453 1.99996580572033 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 14 页，共 29 页 - - - - - - - - - 0.999971141727589 1.9999827300678 迭代次数为 17 由以上输出得到 w取值不同的情况下，得到的满足精度要求的结果，迭代次数分别如输出所示第五次上机作业8. 从函数表x 0.0 0.1 0.2 0.3 0.401 0.5 f(

21、x) 0.39894 0.39695 0.39142 0.38138 0.36812 0.35206 出发，用下列方法计算f(0.15),f(0.31)及 f(0.47)的近似值（1）分段线性插值（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; else fprintf(线性插值 n); for t=1:m z=x

22、(t); 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 15 页，共 29 页 - - - - - - - - - if(zx0(p) fprintf(x%d不在所给自变量范围内，无法进行插值,t); continue; else fori=1:p-1 if(zx0=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(x

23、0,y0,x) 线性插值ans = 0.39404 0.38007 0.35693 即结果为 f(0.15)0.39404,f(0.31) 0.38007,f(0.47) 0.35693 （2）分段二次插值编写的函数如下function R = div2line(x0,y0,x) %分段二次插值p=length(x0); m=length(y0); n=length(x); 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 16 页，共 29 页 - - - - - - - - - i

24、f(p=m) error(输入错误，请重新输入数据); end; for t=1:n if(x(t)x0(p) fprintf(x%d不在所给区间上,t); continue; end; tag=2; m=abs(x(t)-x0(1)+abs(x(t)-x0(2)+abs(x(t)-x0(3); fori=3:p-1 sum=abs(x(t)-x0(i-1)+abs(x(t)-x0(i)+abs(x(t)-x0(i+1); if(sumx0=0.0 0.1 0.195 0.3 0.401 0.5; y0=0.39894 0.39695 0.39142 0.38138 0.36812 0.352

25、06; x=0.15 0.31 0.47; div2line(x0,y0,x) ans = 0.39448 0.38022 0.35725 即分段二次插值方法下，f(0.15) 0.39448,f(0.31) 0.38022,f(0.47) 0.35725 （3）上机编写的程序如下名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 17 页，共 29 页 - - - - - - - - - function R = lagrange(x0,y0,x) %全区间上拉格朗日插值p=lengt

26、h(y0);n=length(x0);m=length(x); %计算函数表和 x 的长度if p = n error(数据输入有误，请重新输入); %若函数表的 x 与y 长度不一致则输入有误else fprintf(拉格朗日插值 n); for t=1:m %利用循环计算每个x 的插值 s=0.0; z=x(t); for k=1:n p=1; fori=1:n ifi=k p=p*(z-x0(i)/(x0(k)-x0(i); end end s=s+y0(k)*p; end %根据拉格朗日插值公式求解y R(t)=s; %输出插值结果end end 上机运行结果为x0=0.0 0.1 0

27、.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; lagrange(x0,y0,x) 拉格朗日插值ans = 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 18 页，共 29 页 - - - - - - - - - 0.39447 0.38022 0.35722 即分段二次插值方法下，f(0.15) 0.39447,f(0.31) 0.38022,f(0

28、.47) 0.35722 9. 解：上机程序如下，为方便起见，将所有操作分在四个函数中进行入口函数function =spline( X,Y,xx,y1_0,y1_18 ) %输出自变量所对应的函数值M=getM(X,Y,y1_0,y1_18);% 先得到 M s=xx; k=length(xx); for a=1:k s(xx(a)=getExp(X,Y,M,xx(a);%计算自变量所在小区间对应曲线的表达式，并根据表达式计算对应的函数值fprintf(s(%d)=%fn,xx(a),s(xx(a); % 输出打印函数值end end 获取 M function M = getM(X,Y,y

29、1_0,y1_1) %得到 M n=length(X); a=0*X;b=a;c=a;h=a;f=a; b=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:n-1); a(n)=1;c(1)=1; Yf(2:n)=Y(2:n)-Y(1:n-1); f(2:n-1)=6*(Yf(3:n)./h(3:n)-Yf(2:n-1)./h(2:n-1)./(h(2:n-1)+h(3:n); f(1)=6*(Yf(2)/h(2)-y1_0)/h(2); f(n)=6*(y1_1-

30、Yf(n)/h(n)/h(n); M=CalM(n,a,b,c,f);% 计算 M end 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 19 页，共 29 页 - - - - - - - - - 计算 M function f = CalM(n,a,b,c,f) % 追赶法求解 M eps=0.1e-15; % 防止参数过小，是的计算误差过大if abs(b(1)eps disp(除数为 0，停止计算 ); return else c(1)=c(1)/b(1); end %追赶法

31、：根据递推算式计算fori=2:n-1 b(i)=b(i)-a(i)*c(i-1); if abs(b(i)X(n5) n1=n5; else n2=n5; end end % % 计算 y y=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); % %结束end 上机试运行，过程如下 X=0.52 3.1 8.0 17.95 28.65 39.62

32、 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.825123 s(4)=10.481311 s(6)=12

33、.363477 s(12)=16.575574 s(16)=19.091594 s(30)=25.386402 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 21 页，共 29 页 - - - - - - - - - s(60)=32.804283 s(110)=36.647878 s(180)=35.917148 s(280)=29.368462 s(400)=16.799784 s(515)=0.542713 根据上述程序运行结果， 可得所有自变量对应的函数值，如上输出结果所示

34、第六次上机作业10. 已知一组实验数据i 1 2 3 4 5 6 7 8 9 ?1 3 4 5 6 7 8 9 10 ?10 5 4 2 1 1 2 3 4 试用最小二乘法求它的多项式拟合曲线，并求出最低点位置。解：试用 matlab 绘图命令，将以上个点绘制在坐标图上，如下图所示名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 22 页，共 29 页 - - - - - - - - - 该函数图像形如二次函数的图像，现使用二次拟合程序如下function xmin,ymin = Se

35、condFitting(x,y) % 最小二乘法，二次拟合 , 形如 a+bx+cx2 x=x; y=y; m=length(x); A(:,1)=ones(m,1); A(:,2)=x; A(:,3)=x.2; cc=Ay; a=cc(1); 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 23 页，共 29 页 - - - - - - - - - b=cc(2); c=cc(3); fprintf(拟合的曲线为 %f%fx+%fx2n,a,b,c); xx=1:0.01:10;

36、yy=cc(1)+cc(2)*xx+cc(3)*xx.2; plot(x,y,r*,xx,yy); ymin=min(yy); c1=a-ymin; p=c b c1; xmin=roots(p); fprintf(最低点的横纵坐标分别为); end 上机运行 x=1 3 4 5 6 7 8 9 10; y=10 5 4 2 1 1 2 3 4; SecondFitting(x,y) 拟合的曲线为 13.459664-3.605309x+0.267571x2 最低点的横纵坐标分别为ans = 6.74 6.7342 图像如下图所示名师资料总结 - - -精品资料欢迎下载 - - - - - -

37、 - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 24 页，共 29 页 - - - - - - - - - 根据程序运行结果，得到拟合方程为y=13.459664-3.605309x+0.267571x2,最低点坐标为（6.74， 6.7342） 。计算方法第七次上机作业15. 用龙贝格算法计算椭圆?2400+?2100 =1 的周长，使误差不超过 10- 4. 解：椭圆周长 L 计算公式为L=4300( ?)2+ 100?20 d ?上机程序如下所示（注：把每次计算结果存储于一个二维数组之中，为了输出如书中所示的表格形式，每一次结果的下标有一定

38、的调整）function = Romberg( fun,a,b,e )%fun 表示被积函数，a 为区间下限，b为上限，e为精度要求%龙贝格算法求椭圆周长名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 25 页，共 29 页 - - - - - - - - - T=zeros();T(1,1)=(b-a)/2*(feval(fun,a)+feval(fun,b);T(2,1)=0.5*T(1,1)+(b-a)/2*(feval(fun,a+(b-a)/2);T(2,2)=(4*T(2

39、,1)-T(1,1)/3;T(3,1)=0.5*T(2,1)+(b-a)/4*(sum(feval(fun,a+(2*(1:2)*(b-a)/4);T(3,2)=(4*T(3,1)-T(2,1)/3;T(3,3)=(16*T(3,2)-T(2,2)/15;T(4,1)=0.5*T(3,1)+(b-a)/8*(sum(feval(fun,a+(2*(1:4)*(b-a)/8);T(4,2)=(4*T(4,1)-T(3,1)/3;T(4,3)=(16*T(4,2)-T(3,2)/15;T(4,4)=(64*T(4,3)-T(3,3)/63;T(5,1)=0.5*T(4,1)+(b-a)/16*(s

40、um(feval(fun,a+(2*(1:8)*(b-a)/16);T(5,2)=(4*T(5,1)-T(4,1)/3;T(5,3)=(16*T(5,2)-T(4,2)/15;T(5,4)=(64*T(5,3)-T(4,3)/63;%以上为求出序列表中前五行的值k=5;while(4*abs(T(k,4)-T(k-1,4)e) k=k+1;T(k,1)=0.5*T(k-1,1)+(b-a)/(2(k-1)*(sum(feval(fun,a+(2*(1:(2(k-2)*(b-a)/(2(k-1);T(k,2)=(4*T(k,1)-T(k-1,1)/3;T(k,3)=(16*T(k,2)-T(k-

41、1,2)/15;T(k,4)=(64*T(k,3)-T(k-1,3)/63;enddisp(4*T);%区间为四分之一圆周，所以乘以4输出为椭圆周长上机运行，结果如下（复制数字结果，排版容易乱，截图如下）运行命令fun=(x)power(300*(sin(x).2)+100),0.5); Romberg(fun,0,pi/2,0.0001) 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 26 页，共 29 页 - - - - - - - - - 由 输 出 结 果 可 得 ， 满

42、足 精 度 要 求 的 椭 圆 的 周 长 计 算 结 果 为96.8845838614657，大约为 96.8846，序列表长宽均为23 计算方法第八次上机作业17.用下列方法求?=23?- 2,? 0,1?0 = 1的数值解（取 h=0.1）,并将计算结果与准确解y= 1 + ?23进行比较 : (1)欧拉方法(2)改进欧拉方法(3)经典 R-K法解：本题 f(x,y)=23?- 2,h=0.1, 上机程序如下function A = Equation(fun,a,b,y0,h)%f表示函数 ,a,b为区间端名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - -

43、- - - - - - - - 名师精心整理 - - - - - - - 第 27 页，共 29 页 - - - - - - - - - %y0为0处的函数值， h为步长% 欧拉方法x=a:h:b;n=(b-a)/h;y=0*x;y(1)=y0;A(1,1)=x(1);A(1,2)=y(1);for k=1:n y(k+1)=y(k)+h*feval(fun,x(k),y(k); A(k+1,1)=x(k+1);A(k+1,2)=y(k+1);end% 改进欧拉方法 A(1,3)=y(1);for k=1:n s1=feval(fun,x(k),y(k);ys=y(k)+h*s1; s2=fe

44、val(fun,x(k+1),ys);y(k+1)=y(k)+0.5*h*(s1+s2);A(k+1,3)=y(k+1);end%R-K方法A(1,4)=y(1);for k=1:n K1=feval(fun,x(k),y(k); K2=feval(fun,x(k)+0.5*h,y(k)+0.5*h*K1); K3=feval(fun,x(k)+0.5*h,y(k)+0.5*h*K2); K4=feval(fun,x(k)+h,y(k)+h*K3);y(k+1)=y(k)+h*(K1+2*K2+2*K3+K4)/6;A(k+1,4)=y(k+1);end% 精确解for k=0:n x=a+k

45、*h;y(k+1)=(1+x2)(1/3);A(k+1,5)=y(k+1);end名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 28 页，共 29 页 - - - - - - - - - end上机运行，结果如下 f=(x,y)2*x*power(y,-2)/3; Equation(f,0,1,1,0.1) 清晰版从左往右，依次是欧拉方法，改进欧拉方法，R-K方法计算的结果名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 29 页，共 29 页 - - - - - - - - -