东南大学数值分析上机作业

习题117.（上机题）舍入误差与有效数（1）编制按从大到小的顺序Sn=、1+，+―+」一，计算»的通用程序。n22-132-1 N2-1(2)编制按从小到大的顺序限=111

H d 1

N2-1(N-l)2-1 22-1,计算％的通用程序。（3）按两种顺序分别计算S］。?，S1q4,S1Q<（4）通过本上机题，你明白了什么？答:1、从大到小顺序：并指出有效位数。（编制程序时用单精度）2、从小到大顺序：1.#include<stdio.h>1.#include<stdio.h>2.intmain()2.intmain()3.{3.{4.floatS=0»TJI»i=2;4.floatS=0,TJI;5.print”，•请输入N：5.print”“请输入II：“）;6.scanf6.7.while(i<=N)7.while(ll>=2)8.{8.{9.T=l/(i*i-l);9.10.i=++i;10.11.S+=T;11.S+=T;12.)12.}13.printf(MS=%f\n'\S);13.printf(MS=%f\n”>S);14.system("pause");14.system("pause");15.return0;15.return0;16.}16.}3、计算结果:N102104106从大到小Sn0.7400490.7498520.749852从小到大Sn0.7400500.7499000.749999精确值0.7400490.74990.749999有效位数从大到小Sn643从小到大Sr5464、感想：算法的不同导致计算结果的误差不同，一定要好好设计算法，尽量使算法精确。通过上面的例子，可以看出按不同的计算顺序所得到的结果是不同的，按从大到小的顺序计算的值

与精确值有较大的误差，而按从小到大的顺序计算的值与精确值吻合。计算机在进行数值计算时会出现“大数吃小数”的现象，导致计算结果的精度有所降低，使用计算机进行同号数的加法时，采用绝对值较小者先加的算法，得到的结果的相对误差较小。5、附图：0ll.cttinclude<stdio.h>int0ll.cttinclude<stdio.h>intmain()■C:\Windows\$y$tem32\Debug\ll.exe请输入N：100S=0.740049请按任意键继续..FloatS=0,T,N»i=2;print请输入N：scanf(”F1&N);while(i<=N)S*=T;printf(,eS-%F\nie,S);systemC^pause**);return0;图1从大到小算法及结果■C:\Windows\system32\Debug\ll.exel青输入■C:\Windows\system32\Debug\ll.exel青输入N：1。。S-0.740050请按任意键继续.Itinclude<stdio.h>intmain(){FloatS=0,T,N;printfL请榆入N：scanF(e,V»&M);while(N>-2)<T-1/(H«N-1);N H;S*-T；>printFC'S-%F\n"SS);systemC^ause**);return0;图2从小到大算法及结果习题220.（上机题）Newton迭代法（1）给定初值4及容许误差£,编制Newton法解方程f（x）=0根的通用程序。（2）给定方程f（x）=（/3-乂=0,易知其有三个根xf ,xj=0,K=W0

.由Newton方法的局部收敛性可知存在b>0,当％£（-5,磔时，Newton迭代序列收敛于根亡。试确定尽可能大的.试取若干初始值，观察当圣£（一8,-1）,（-1,-J）,（―色力，（3,1）,（1,8）时Newton序列是否收敛以及收敛于哪一个根。（3）通过本上机题，你明白了什么？答:1、通用程序：.#include<stdio.h>.1、通用程序：.#include<stdio.h>.#include<math.h>p.floatf(floatx){floatf;f=x*x*x/3-x;return(f);}floatdf(floatx)1floatdf;df=x*x-l;return(df);}115.intmain()T6.{intk=0;printfC,请输入初值：”)；scanf("%fM,&x0);TOC\o"1-5"\h\zdo{d=-f(x0)/df(x0);xl=x0*-d;k++;x0=xl;}while(fabs(d)>eps);print",迭代％d次后，根=%八n”，k,x0);system("pause");return0;}#include<stdio.h>.#include<math.h>floatf(floatx)floatf;f=x*x*x/3-x;return(f#include<stdio.h>.#include<math.h>floatf(floatx)floatf;f=x*x*x/3-x;return(f);}float}floatdf(floatx)1floatdf;df=x*x-l;return(df);1115.intmain(){floatx0=0,xl>d;floateps=0.000005,delta=0;while(x0==0){delta=delta+0.000001;. x0=x0^delta;do1d=-f(x0)/df(x0);xl=x&i-d;x0=xl;}. while(fabs(d)>eps);}delta=delta-0.000001;printf(“尽可能大的值=%八rTjdelta);. system。'pause");return0;}运算结果为：0.7745963、测试是否在给定区间收敛：在五个区间内分别取初值-100000,-0.8,0.5,0.9,110000测试是否收敛及收敛到哪个根，结果如下表：初值-100000-110000是否收敛是是是是是收敛于-V3-V30V3V34、感想：选取的迭代初值不同，则迭代序列将有可能收敛于不同的根，牛顿法对迭代初值的要求较高。有些区间收敛范围相当大，有些却非常小。另外，我发现容许误差的大小对迭代结果的精度有很大的影响，如果想得到比较精确的结果，要将容许误差设置小一点。5、附图:目lx lol-—return(F);)floatdf(floatx)Floatdf;df-xKX-1;return(df);)dntmain()(floatx0,x1fdteps=0.0000005;intk-A;print”“话输入初值一,)；scanF(,V\fo(0);do<d-f(x0)/df(x0>;x1=x0*d;…；xb-xl;>while(fabs(d)>eps);pTn"(,・迭代td次乳根T『\n・,,3町;system"pause");return0;■C:\w\do*w5\5y$tem32\Debug\Lexe79?LT于旧土.,|江口照।9；后■根:1.732。51幽拄续・・・・।C:\vjindows\3ystem32\Debug\l.exe।C:\vjindows\3ystem32\Debug\l.exe图1利用牛顿法求方程的近似解目1.CFloatdf;df=x*x-1;return(dF);}intmain()<Floatx0=8»x1,d;Floateps=0.O00Q05,delta=O;while(x0==9){delta=delta*8.800801;x0=xO+delta;do<d=-f(xO)/dF(xfl);x1=x8+d;X0=x1；>while(Fabs(d)>pps);>delta=delta-8.000001;printer尽可能大的值咤F\n“,delta);systemCpause");return0;图2确定收敛范围习题339.(上机题)列主元Gauss消去法对于某电路的分析，归结为求解线性方程组RI=V，其中R为9阶矩阵，参见课本127页Vr=(-15,27,-23,0,-20,12,-7,7,0)(1)编制解n阶线性方程组Ax=b的列主元三角分解法的通用程序；(2)用所编制的程序解线性方程组RI=V,并打印出解向量，保留五位有效数:(3)本编程之中，你提高了哪些编程能力？筝1、通用程序■#include<stdio.h>#include<math.h>p.#definedeltale-64.#defineN100main(){int j〉t〉r“n»u»c=0;float p,L»max,s;float X[N];float a[N][N+l];printf(“请输入系数矩阵的阶数：\n”)；scanf("%d",&n);printf(“输入的增广矩阵系数，中间用空格隔开：\n”)；for(i=0;i<n;i++)for(j=0;j<n+l;j++)scanf(M%fM,&a[i][j]);printfC,你输入的增广矩阵为：\n”)；for(i=0;i<n;i++)TOC\o"1-5"\h\z{for(j=0;j<n+l;j++)print"八％八printf}for(j=0;j<n-l;j++){{max=fabs(a[j][j]);r=j;}for(i=j+l;i<n;i++)if(fabs(a[i][j])xnax){max=a[i][j];r=i;・・}if(fabs(a[i][j])<delta)printf("矩阵奇异,,)；for(t=j;t<n+l;t++)|P=a[j][t];a[j][t]=a[r][t];a[r][t]=p;}for(i=j+l;i<N;i++)|L-a[i][j]/a[j][j];for(t=j;t<n+l;t++)a[i][t]=a[i][t]-L*a[j][t]；}}printf("输出原方程的解为：\n”)；X[n-l]=a[n-1][n]/a[n-1][n-1];for(i=n-2;i>=0;i--){s=a[i][n];for(j=i+l;j<n;j++)s=s-a[i][j]*X[j];X[i]=s/a[i][i];}for(u=0;u<n;u++)for(j=0;j<n+l;j++){printf("%12f",a[u][j]);C++；if(c%(rH-l)==0)printf("\n");}for(i=0;i<n;i++)|printf(',x(%d)=%.4f\n->i+l>X[i]);if(i==n-1)printf("An'");}system("pause");return0;2、运行程序，方程的解为:

xlx2x3x4x5x6x7x8x9-0.28923034544-0.71281-0.22061-0.430400.15431-00578230.201050.290233、感想：首先，通过本上机题，提高了我运用循环语句的能力；其次，是我对各种排序算法有了更深的理解：最后，是我对二维数组有了更好的把握，并且认识到MATLAB在处理矩阵运算上的优势。在上机的过程中，加深了对Gauss消去法的认识。4、附图：C:\vjindows\system32\Debug\tes.exeC:\vjindows\system32\Debug\tes.exe请输入系数矩阵的阶数:A000-304100你输入的增广矩阵为:31.00-13.000.000-50027-270W27W-931TSWWWWW-23输入的增广矩阵系数，中间用空格隔开：31-13WWW-1MWWW-15-1335-9W-11WW0-10.000.00,UMU.UUU.MM-9.MWM.WW23.000.000.00-IM.MM7V.SS000.000.000.00W.MMW.MW0.00-5.000.0020.007.00-2.000.0027.00-2.0029.0019.000.000.00-5.00A000-304100你输入的增广矩阵为:31.00-13.000.000-50027-270W27W-931TSWWWWW-23输入的增广矩阵系数，中间用空格隔开：31-13WWW-1MWWW-15-1335-9W-11WW0-10.000.00,UMU.UUU.MM-9.MWM.WW23.000.000.00-IM.MM7V.SS000.000.000.00W.MMW.MW0.00-5.000.0020.007.00-2.000.0027.00-2.0029.0019.000.000.00-5.00-9.00运算后的矩阵及方程的解为:31.000000 -13.0000009.00.909000000-10,0000000.000000 0.000000 0.000000-15.000000A.AA29.548388A.AAUSSSU.UUUUUUU.UUUUUU 2U.7096770.0000000.000000-10.0000000.000.0000000.000000-10.0000000.009000 0.000000 0.000000 -16.6921410.000000 0.000000 -0.000000 75.4G1273 -31.185629 -0.451999 0.000000 0.000000 -9.000000 -5.906896图1输入增广矩阵printer运算后的矩阵及方程的解为：\n")printer运算后的矩阵及方程的解为：\n");X[n-1]-a[n-1][n]/a[n-1][n-1];for(i=n-2;i>=fl;i-)<s=a[i][n];for(J=1*1;s=5-a[i][j]*X[j];X[i]-5/a[i][i];)for(u=n；u<n;u+4)€or(j=9;j<n*1;j++)<printfC%12F-,a[u][j]);iF(c%(n*1)==0)printfCAn1');>for(i=0;l<n;i**)<printF(ilx(^d)-%.5F\n'i,i*1,X[i]);if(i==n-1)printf(M\n");)systemC'pausereturn0;，C\windows\system32\Dcbug\tcs.exe0.000000 0.000000 28.258734-10J0000 0.000000 0.000000 -16.6921410.WWWWHW M.MWMWMH -H.MMMWWW 75HUWS W.H0OOOO -9.WUU00O -5・，96B960•000009e.000900 0.000000 -0.(9000 -5.000000 -3.577994 -22.3483160.000000 0.000800 0.000000 8JGQQG -0.784718 -0.561543 8.492S740.WWWWHW M.MWMWMH H.MMMWWW WJB697 -M.513187 -0.367236 -1.4469540•000009 e.000900 0.000000 0.19880 26.413086 2.417996 4.6081010.000000 0.000800 0.000000 -0.(0000 0.000000 27.389503 7.949220x<l>="0.28923x<2)=0.34544x<3>0.71281x<4>-0.22061x<5)0.43010x<6)=0.1S431x<?)=-0.05782x<«)=M.2WlkJ5x<95-0.29023清按任意键维续・•• 习题437.(上机题)3次样条插值函数(1)编制求第一型3次样条插值函数的通用程序;(2)已知汽车门曲线型值点的数据加下：1012345678910X1012345678910Yi2.513.304.044.705.225.545.785.405.575.705.80端点条件为y；=0.8,y；o=O.2,用所编程序求车门的3次样条插值函数S(x),并打印出S(i+0.5),i=0,1,…，9o答:(1)通用程序:#include<stdio.h>#include<math.h>#include<stdlib.h>#definedeltale-6#defineN100intmain(){intfloatp,L»max》s,h»Dl»D2;floatx[N],M[N],y[H];floata[N][N+l]={0};/嘛入插值条I牛*/printf(“请输入插值节点的个数：\n”)；scanf(M%d'\&n);. printf(“请从左到右输入插值节点：\n");for(i=0;i<n;i++)TOC\o"1-5"\h\z1scanf(M%f'\&x[i]);}printf(“请输入插值节点对应的函数值：\n");for(i=0;i<n;i++){scanf(M%f'\&y[i]);}printf(“请输入插值初始条件：\n“)；printf(My,(0)=M)；scanf(M%f'\&Dl);printf(My,(n)=M)；scanf(M%f'\&D2);/*产生弯矩方程组增广施阵系数*/a[0][l]=l；a[n-1][n-2]=l;for(i=0;i<n;i++)a[i][i]=2;for(i=l;i<n-l;i++)TOC\o"1-5"\h\z{h=(x[i]-x[i-l])/(x[i+l]-x[i-l]);a[i][i-l]=h;a[i][i+l]=l-h;}a[0][n]=6*((y[n-2]-y[n-l])/(x[n-2]-x[n-l])-D2)/(x[n-2]-x[n-l]);a[n-l][n]=6*((y[l]-y[0])/(x[l]-x[0])-Dl)/(x[l]-x[0]);for(i=l;i<n-l;i++){. a[i][n]=6*((y[i+l]-y[i])/(x[i+l]-x[i])-(y[i]-y[i-l])/(x[i]-x[i-l]))/(x[l+l].x[i.l]);}print,请确认弯矩方程的增广矩阵为：\n“);for(i=0;i<n;i++){for(j=0;j<n+l;j++)printf(“％3.2f\t，a[i][j]);printf}for(j=0;j<n-l;j++){{max=fabs(a[j][j]);r=j;}for(i=j+l;i<n;i++)if(fabs(a[i][j])xnax){max=a[i][j];r=i;}if(max<delta)printfC,矩阵奇异\n“)；for(t=j;t<n+l;t++){P=a[j][t];a[j][t]=a[r][t];a[r][t]=p;}for(i=j+l;i<N;i++)TOC\o"1-5"\h\z{L=a[i][j]/a[j][j];for(t=j;t<n+l;t++)a[i][t]=a[i][t]-L*a[j][t];}}/嗝出弯矩方程的计算结果*/printf(••整理后的增广矩阵为：\n”)；M[n-l]=a[n-1][n]/a[n-1][n-1];for(i=n-2;i>=0;i--)1s=a[i][n];. for(j=i+l;j<n;j++)s=s-a[i][j]*M[j];M[i]=s/a[i][i];TOC\o"1-5"\h\z}for(u=0;u<n;u++)for(j=0;j<n+l;j++)1printf(M%3.2f\t'\a[u][j]);. C++;if(c%(n+l)==0)printf(M\n");}printf(“弯矩方程解向量为：\n”)；for(i=0;i<n;i++){printf(MM(%d)=%f\n"/i+lJ4[i]);if(i==n-l)printf(M\nM);}/*计算插值函数*/printf(“样条函数如下：\n”)；for(i=0;i<n-l;i++){printf(MS(x)=%f+%f(x-%.2f)+%f(x-%.2f)A2+%f(x-%.2f)A3\t%2.1f<x<%2.1f，，,y[i],((y[i+l]-y[i])/(x[i+l]-x[i])-(M[i]/3+M[i+l]/6)*(x[i+l]-x[i])),x[i],M[i]/2,x[i],(M[i+l]-M[i])/6/(x[i+l]-x[i]),x[i],x[i],x[i+l]);printf(M\n");}printf(，3(1+0.5)的值如下：\n”)；for(i=0;i<n-l;i++){

printf(MS(%d+0.5)=%f\iJy[i]+((y[i+l]-y[i])/(x[i+l]-x[i])-(M[i]/^-M[i+l]/6)*(x[i+l]-x[i]))*(i+0.5-x[i])+M[i]/2*(i+0.5-x[i])*(i+0.5-x[i])+(M[i+l]-M[i])/6/(x[i+l]-x[i])♦(!+©.5-x[i])*(i+0.5-x[i])*(i+0.5-x[i]));printf("\n");}. system(''pause");return0;|121.)(2)打印出3次样条插值函数S(x),并打印出S(i+0,5),(i=0,1,…，9)1、样条函数如下：' 2.510000+0.690000X+0.189038X2-0.089039X3 (0,1)S(x)=(3)附图:3.300000+0.800961(x-1)-0.078078(x-I)2+0.017117(x-I)3(1,2)4.040000+0.696156(x-2)-0.026728(x-2)2-0.009428(x-2)3(2,3)4.700000+0.614416(x-3)-0.055011(x-3)2-0.039405(x-3)3(3,4)5.220000+0.386178(x-4)-0.173227(x-4)2+0.107049(x-4)3(4,5)5.540000+0.360870(x-5)+0.147919(x-5)2-0.268789(x-5)3(5,6)5.780000-0.149659(x-6)-0.658448(x-6)2+0.428107(x-6)3(6,7)5.400000-0.182234(x-7)+0.625873(x-7)2-0.273639(x-7)3(7,8)5.570000+0.248596(x-8)-0.195043(x-8)2+0.076447(x-8)3(8,9)<5.70000+0.08785l(x-9)+0.034299(x-9)2-0.022149(x-9)S(x)=(3)附图:X.5Yi2.89113.68314.38024.98855.38325.72385.59415.43115.65515.74972、S(i+0.5)的值如下:■*XAUscrsSzh^nah3n\Oo<imcnts\CFrccXProjects,vi11请从左利右脸人科侦盲点:0访箱入插值节点对应的由数值,2.513.3I.011.75.225.$15・785.57喻输入插值切始条件；/⑹-0.87,(n)=0.2中文 QQ府自购入次T:为:图1 输入插值条件

0.0.弯的)=)=认后^^^ooooooooooow^ooooooooooo((5050000000000-0000000000082矩L.弯的)=)=认后^^^ooooooooooow^ooooooooooo((5050000000000-0000000000082矩L2.0.0.s0.().().0.0.0.isLL-0O.O.O.O.O.O.O.O.方0000500000000000000000广007500000000000000000增..0.0.0.为O.SLYO.O.O.O.O.O.S广0050005000000000000000为..0.0.0.O00(50005000000000000Oooooooooooo0005050000.S0.0.0.O00600000000058•ooooooO••••O ooO1oooooOO000^/00000000058•oooooO00007000000000580000000.0.SO.LO.0.O.O.O.O.0.000.000.000.000.000.000.600.000.000.000.000.000.0()-0.150.000.000.000.000.000.00-0.240.000.000.000.000.000.00-0.420.500.000.000.000.000.00-0.602.000.500.000.000.000.00-0.240.502.000.500.000.000.00-1.860.000.502.000.500.000.001.650.000.000.502.000.500.00-0.120.000.000.000.502.000.50-0.090.000.000.000.001.002.00-0.060.000.000.000.000.000.000.600.000.000.000.000.000.00-0.300.000.000.000.000.000.00-0.150.000.000.000.000.000.00-0.380.500.000.000.000.000.00-0.501.870.500.000.000.000.00-0.11-0.001.870.500.000.000.00-1.830.00-0.001.870.500.000.002.140.000.00-0.001.870.500.00-0.690.000.000.00-0.001.870.500.100.000.000.000.00-0.001.73一0.11图2 生成弯矩方程的增广矩阵d(7)=-l.3168964(8)=1.251746M⑼二-0.390087W(10)=0.068597M(U)=-0.064299样条函数如下：0.0<x<1.01.0<x<2.02.0<x<3.03.0<x<4.00<x<5.00<x<6.00<x<7.00<x<8.00<x<9.00<x<10.0S(x)=2,510000+0.0.0<x<1.01.0<x<2.02.0<x<3.03.0<x<4.00<x<5.00<x<6.00<x<7.00<x<8.00<x<9.00<x<10.0S(x)=3.300000+0.800961(x-1.00)f078078(xT.00厂2+0.017U7(x-l.00厂3S(x)=4.700000+0.614416(x-3.00)S(x)=4.700000+0.614416(x-3.00)+-0.05501l(x-3.00厂2+-0.039405(x-3.00)*3S(x)=5.220000+0.386178(x-4.00)+-0.173227(x-4.00)*2+0.107049(x-4.00”3S(x)=5.540000+0.360870(x-5.00)+0.147919(x-5.00)'2+-0.268789(x-5.00)’3>(x)=5.780000+-0.149659(x-6.00)+-0.658448(x-6.0042+0.428107(x-6.00厂35(x)=5.400000十-0.182234(x-7.00)+0.625873(x-7.00厂2十-0.273639(x-7.00厂35(x)=5.570000+0.248596(x-8.00)+-0.195043(x-8.00厂2+0.076447(x-8.00厂3S(x)=5.700000+0.087851(x-9.00)+0.034299(x-9.00)A2+-0.022149(x-9.00)’35(i+0.5)的值如下:S(0+0.5)=2.891130S(l+0.5)=3.6831015(240.5)=4.3802175(3+0.5)=4.988530式4+0.5)=5.383163式5+0.5)=5.7238165(6+0.5)=5.594072S(7+0.5)=5.431146S(8+0.5)=5.655093S(9+0.5)=5.749731图3 输出插值函数及计算结果习题623.（上机题）常微分方程初值问题数值解常微分方程初值问题数值解（1）编制RK4方法的通用程序；（2）编制AB4方法的通用程序（由RK4提供初值）；（3）编制AB4-AM4预测校正方法通用程序（由RK4提供初值）；（4）编制带改进的AB4-AM4预测校正方法通用程序（由RK4提供初值）；（5）对于初值问题g=-%2y2Iy（o）=3取步长h=01,应用（1）〜（4）中的四种方法进行计算，并将计算结果和精确解y（x）=3/（l+%3）作比较;（6）通过本上机题，你能得到哪些结论？答:RK4通用程序:#include<stdio.h>#include<stdlib.h>#include<math.h>#defineH0.1intmain(){floatx»y,r=3;floatkl,k2,k3,k4;floatfuni(floatyfloat);inti»n;printfC,请输入x,y的初值：\n“)；scanf(M%f%f\&x,&y);printf(・.请输入运算次数：\n“)；scanfC'Xd'\&n);printf(“运算结果为：\n”)；for(i=0;i<n;i++){printf(",x%-2d=%.2f>\ty%-2d=%f,\ty(%-2d)=%fJ\tr(%-2d)=%f'\i.x>i>y>iJrprintf("\rT);kl=funl(x>y);k2=funl(x+H/2,y+H*kl/2);k3=funl(x+H/2,y+H*k2/2);k4=funl(x+H,y+H*k3);x=x+H;y=y+H*(kl+2*k2+2*k^k4)/6;r=3/(l+x*x*x);1system(Mpause");return0;}float funl(float x^floaty){floatf;f=-x*x*y*y;returnf;}AB4通用程序：#include<stdio.h>.#include<stdlib.h>3・#include<math.h>#defineH0.1intmain(){floatx»y〉r=3;floatkl,k2,k3>k4>t[100]={0};floatfuni(float,float);intprintfC,请输入x,y的初值：\n“)；scanf(M%f%f\&x,&y);print£(,,请输入运算次数：\n“)；scanfC^d'\&n);printf(”运算结果为：\nM);for(i=4;i<n;i++)TOC\o"1-5"\h\z{if(i<5){for(j=0;j<4;j++){printf(Mx%-2d=%.2f,\ty%-2d=%f,\ty(%-2d)=%fMJ.xJ.yJ.r);printf("\n");kl=funl(x,y);k2=funl(x+H/2,y+H*kl/2);k3=funl(x+H/2,y+H*k2/2);k4=funl(x+H»y+H*k3);t[j]=y;x=x+H;y=y+H♦(kl+2*k2+2*k3+k4)/6;r=3/(l+x*x*x);}x=x-H;}t[i]=t[i-l]+H*(55Hunl(xj[i-l])-59*funl(x-H,t[i-2])+37Hunl(x-2*H,t[i-3])-9*funl(x-3*H,t[i-4]))/24;x=x+H;r=3/(l+x*x*x);printf(',x%-2d=%.2f>\ty%-2d=%f,\ty(%-2d)=%f,\tr(%-2d)=%fM>i>x>i>t[i]>i,r,i,r-t[i]);printf(”\n”)；}system(Mpause");return0;}float funl(float x^floaty){floatf;f=-x*x*y*y;returnf;}AB4-AM4通用程序：#include<stdio.h>.#include<stdlib.h>#include<math.h>#defineH0.1intmain(){float x»y,r=3;float kl,k2,k3>k4>t[100]={0}>Y[100]={0};float funi(float,float);intprintf(“请输入x,y的初值：\n“)；scanf("%f%f”>&x〉&y);printf(“请输入运算次数：\n");scanfprintf(•,运算结果为：\n“);for(i=4;i<n;i++)TOC\o"1-5"\h\z{if(i<5){for(j=0;j<4;j++)1printf("x%-2d=%.23\ty%・2d=%f,\ty(%-2d)=%f"，j,x,j,y,j,r);第16页printf("\n");kl=funl(x>y);k2=funl(x+H/2,y+H*kl/2);k3=funl(x+H/2,y+H*k2/2);k4=funl(x+H»y+H*k3);Y[j]=y;x=x+H;y=y+H*(kl+2*k2+2*k3+k4)/6;r=3/(l+x*x*x);}x=x-H;}t[i]=Y[i-l]+H*(55*funl(x>Y[i-l])-59*funl(x-H,Y[i-2])+37*funl(x-2*H>Y[i-3])-9*funl(x-3*H,Y[i-4]))/24;Y[i]=Y[i-l]+H*(9*funl(x+H,t[i])+19*funl(x,Y[i-l])-5*funl(x-H,Y[i-2])+funl(x-2*H,Y[i-3]))/24;. x=x+H;r=3/(l+x*x*x);printf(',x%-2d=%.2f>\ty%-2d=%f,\ty(%-2d)=%fJ\tr(%-2d)=%f'\i.x>i,Y[i]>i,r,i,r-Y[i]);printf("\rT);}. system(Mpause");return0;}. float funl(float x^floaty){. floatf;f=-x*x*y*y;returnf;}(4)改进的AB4-AM4通用程序：#include<stdio.h>.#include<stdlib.h>#include<math.h>#defineH0.1intmain(){|7. float x»y,r=3;float kl,k2Jk3,k4,t[100]={0}>Y[100]={0}JR[100]={0};float funi(float,float);intprintfC,请输入x,y的初值：\n“)；scanf(”%f%f”,&x,&y);printf(,,请输入运算次数：\n“)；scanfC'%d'\&n);printf(••运算结果为：\n");for(i=4;i<n;i++)TOC\o"1-5"\h\z{if(i<5){for(j=0;j<4;j++){printf("x%-2d=%.2f,\ty%-2d=%f,\ty(%-2d)=%fJ,xJ.yJ,r);printf(M\n");kl=funl(x,y);k2=funl(x+H/2,y+H*kl/2);k3=funl(x+H/2,y+H*k2/2);k4=funl(x+H>y+H*k3);R[j]=y;x=x+H;y=y+H*(kl+2*k2+2*k3+k4)/6;r=3/(l+x*x*x);}x=x-H;}t[i]=R[i-l]+H*(55*funl(x,R[i-l])-59*funl(x-H,R[i-2])+37*funl(x-2+H,R[i-3])-9*funl(x-3*H,R[i-4]))/24;Y[i]=R[i-l]+H*(9*funl(x+H,t[i])+19*funl(x,R[i-l])-5*funl(x-H,R[i-2])+funl(x-2*H,R[i-3]))/24;R[i]=251.0/270*Y[i]+19.0/270*t[i];x=x+H;r=3/(l+x*x*x);printf(',x%-2d=%.2f,\ty%-2d=%f,\ty(%-2d)=%f,\tr(%-2d)=%f'\i.x>i>R[i]>printf("Xn");}system(Mpause");return0;1float funl(float x^floaty){floatf;f=-x*x*y*y;returnf;51・ }

(5)数据比较：•C:\Users\zhanghan\Documents\C-Free\P9jects\3ciYTCZ\RK4.exe'请输入运算次数:116运算结果为：xO=0.请输入运算次数:116运算结果为：xO=0.00,xl=0.10,x2=0.20,x3=0.30,x4=0.40,二0.50,x6=0.60,x7=0.70,x8=0.80,x9=0.90,x10=1.00,xll=l.10,xl2=1.20,xl3=1.30,xl4=1.40,x15=1.50,请按任意键继续.yO=3.000000,yl=2.997003,y2=2.976190,y3=2.921129,y4=2.818389,y5=2.664672,y6=2.465203,y7=2.233079,y8=1.984950,y9=1.737043,yl0=l.502194,y11=1.288763,yl2=l.100724,yl3=0.938710,yl4=0.801135,yl5=0.685334,y(0)=3.000000y(l)=2.997003y(2)=2.976191y(3)=2.921129y(4)=2.819549,y(5)=2.666667,y(6)=2.467105,y(7)=2.233805,y(8)=1.984127,y(9)=L735107,y(10)=L500000,y(ll)=L28700Ly(12)=1.099706,y(13)=0.938379,y(14)=0-801282,y(15)=0.685714,r(4)=0.001160r(5)=0.001994r(6)=0.001903r(7)=0.000726r(8)-0.000824r(9)=-0.001936r(10)=-0.002195r(ll)=-0.001762r(12)=-0.001017r(13)=-0.000331r(14)=0.000147r(15)=0.000380请输入x,y的初值:03请输入运算次数：16运算结果为：xO=0.00,yO=3.000000,y(0)=3.000000,r(0)=0.000000xl=0.10,yl=2.997003,y(l)=2.997003,r(l)=0.000000x2=0.20,y2=2.976190,y(2)=2.976191,r(2)=0.000000x3=0.30,y3=2.921129,y(3)=2.921129,r(3)=0.000001x4=0.40,y4=2.819547,y(4)=2.819549,r(4)=0.000002x5=0.50,y5=2.666663,y(5)=2.666667,r(5)=0.000003x6=0.60,y6=2.467100,y(6)=2.467105,r(6)=0.000005x7=0.70,y7=2.233799,y(7)=2.233805,r(7)=0.000006x8=0.80,y8=1.984123,y(8)=1.984127,r(8)=0.000004x9=0.90,y9=L735107,y(9)=1.735107,r(9)=-0.000000xl0=l.00,yl0=l.500006,y(10)=1.500000,r(10)=-0.000006xll=l.10,yl1=1.287012,y(11)=1.287001,r(ll)=-0.000011xl2=1.20,yl2=l.099722,y(12)=1.099706,r(12)=-0.000016xl3=1.30,yl3=0.938397,y(13)=0.938379,r(13)=-0.000018xl4=l.40,yl4=0.801300,y(14)=0.801282,r(14)=-0.000018xl5=1.50,yl5=0.685732,y(15)=0-685714,r(15)=-0.000018请按任意键继续.••图1RK4解法产'•C:\Users\zhanghan\Documents\C-Free\Projects\3ciYTCZ\AB4.exe，请输入X,y的初值:0图2AB4解法zT"C:\Users\zhanghan\Document