东南大学数值分析上机实验题下.docx

1、东南大学数值分析上机实验题下 数值分析上机报告 姓名： 学号： 2013年12月22日第四章38.（上机题）3次样条插值函数（1）编制求第一型3次样条插值函数的通用程序; （2） 已知汽车曲线型值点的数据如下：0123456789102.513.304.044.705.225.545.785.405.575.705.80端点条件为=0.8， =0.2。用所编制程序求车门的3次样条插值函数S(x),并打印出S(i+0.5)(i=0,1,9)。解：（1）#include #include float x1100,f1100,f299,f398,m100,a100101,x,d100;float c

2、99,e99,h99,u99,w99,y_0,y_n,arr,s;int i,j,k,n,q;void selectprint(float y) if (y0)&(y!=1) cout+y; else if (y=1) cout+; else if (y0) couty; void printY(float y) if (y!=0) couty; float calculation(float x) for (j=1;j=n;j+) if (x=x1j) s=(float)(f1j-1+cj-1*(x-x1j-1)+mj-1/2.0*(x-x1j-1)*(x-x1j-1)+ej-1*(x-x1

3、j-1)*(x-x1j-1)*(x-x1j-1); break; return s;void main() do coutn; if (n100)|(n1) cout请重新输入整数(1.100):100)|(n1); cout请输入Xi(i=0,1,.,n):; for (i=0;ix1i; coutendl; cout请输入Yi(i=0,1,.,nn):; for (i=0;if1i; coutendl; couty_0y_n; coutendl; for (i=0;in;i+) hi=x1i+1-x1i; for (i=1;in;i+) ui=(float) (hi-1/(hi-1+hi)

4、; for (i=1;in;i+) wi=(float) (1.0-ui); for (i=0;in;i+) f2i=(f1i+1-f1i)/hi; /一阶差商 for (i=0;in-1;i+) f3i=(f2i+1-f2i)/(x1i+2-x1i); /二阶差商 for (i=1;in;i+) di=6*f3i-1; /求出所有的di d0=6*(f20-y_0)/h0; dn=6*(y_n-f2n-1)/hn-1; for (i=0;i=n;i+) for (j=0;j=n;j+) if (i=j) aij=2; else aij=0; a01=1; ann-1=1; for (i=1;

5、in;i+) aii-1=ui; aii+1=wi; for (i=0;i=n;i+) ain+1=di; for (i=1;i=n;i+) /用追赶法解方程,得Mi arr=aii-1; for (j=0;j=0;i-) mi=(ain+1-aii+1*mi+1)/aii; for (i=0;in;i+) /计算S（x）的表达式 ci=(float) (f2i-(1.0/3.0*mi+1.0/6.0*mi+1)*hi); for (i=0;in;i+) ei=(mi+1-mi)/(6*hi); for (i=0;in;i+) coutX属于区间x1i,x1i+1时endlendl; cout

6、S(x)=; printY(f1i); if (ci!=0) selectprint(ci); cout(x; printY(-x1i); cout); if (mi!=0) selectprint(float)(mi/2.0); for (q=0;q2;q+) cout(x; printY(-x1i); cout); if (ei!=0) selectprint(ei); for (q=0;q3;q+) cout(x; printY(-x1i); cout); coutendlendl; cout针对本题，计算S(i+0.5)(i=0,1,.,9):endl;for (i=0;i10;i+)

7、 if (i+0.5=x10) calculation(float)(i+0.5); coutS(i+0.5)=sendl; else couti+0.5超出定义域endl; ;coutendl;（2）编制的程序求车门的3次样条插值函数S(x)：x属于区间0,1 时；S(x)=2.51+0.8(x)-0.0014861(x)(x)-0.00851395(x)(x)(x)x属于区间1,2 时；S(x)=3.3+0.771486(x-1)-0.027028(x-1)(x-1)-0.00445799(x-1)(x-1)(x-1)x属于区间2,3 时；S(x)=4.04+0.704056(x-2)-0

8、.0404019(x-2)(x-2)-0.0036543(x-2)(x-2)(x-2)x属于区间3,4 时；S(x)=4.7+0.612289(x-3)-0.0513648(x-3)(x-3)-0.0409245(x-3)(x-3)(x-3)x属于区间4,5 时；S(x)=5.22+0.386786(x-4)-0.174138(x-4)(x-4)+0.107352(x-4)(x-4)(x-4)x属于区间5,6 时；S(x)=5.54+0.360567(x-5)+0.147919(x-5)(x-5)-0.268485(x-5)(x-5)(x-5)x属于区间6,7 时；S(x)=5.78-0.14

9、9051(x-6)-0.657537(x-6)(x-6)+0.426588(x-6)(x-6)(x-6)x属于区间7,8 时；S(x)=5.4-0.184361(x-7)+0.622227(x-7)(x-7)-0.267865(x-7)(x-7)(x-7)x属于区间8,9 时；S(x)=5.57+0.256496(x-8)-0.181369(x-8)(x-8)+0.0548728(x-8)(x-8)(x-8)x属于区间9,10 时；S(x)=5.7+0.058376(x-9)-0.0167508(x-9)(x-9)+0.0583752(x-9)(x-9)(x-9)S(0.5)=2.90856

10、S(1.5)=3.67843 S (2.5)=4.38147S(3.5)=4.98819 S(4.5)=5.38328 S(5.5)=5.7237S(6.5)=5.59441 S(7.5)=5.42989 S(8.5)=5.65976S(9.5)=5.7323第六章21（上机题）常微分方程初值问题数值解（1）编制RK4方法的通用程序；（2）编制AB4方法的通用程序（由RK4提供初值）；（3）编制AB4-AM4预测校正方法的通用程序（由RK4提供初值）；（4）编制带改进的AB4-AM4预测校正方法的通用程序（由RK4提供初值）;（5）对于初值问题取步长,应用（1）（4）中的四种方法进行计算，并将

11、计算结果和精确解作比较；（6）通过本上机题，你能得到哪些结论？解：#include#include#include#includeofstream outfile(data.txt);/此处定义函数f(x,y)的表达式/用户可以自己设定所需要求得函数表达式double f1(double x,double y) double f1; f1=(-1)*x*x*y*y; return f1;/此处定义求函数精确解的函数表达式double f2(double x) double f2; f2=3/(1+x*x*x); return f2;/此处为精确求函数解的通用程序void accurate(do

12、uble a,double b,double h) double x100,accurate100; x0=a; int i=0; outfile输出函数准确值的程序结果:n; do xi=x0+i*h; accuratei=f2(xi); outfileaccuratei=accuratein; i+; while(i(b-a)/h+1);/此处为经典Runge-Kutta公式的通用程序void RK4(double a,double b,double h,double c) int i=0; double k1,k2,k3,k4; double x100,y100; y0=c; x0=a;

13、 outfile输出经典Runge-Kutta公式的程序结果:n; do xi=x0+i*h; k1=f1(xi,yi); k2=f1(xi+h/2),(yi+h*k1/2); k3=f1(xi+h/2),(yi+h*k2/2); k4=f1(xi+h),(yi+h*k3); yi+1=yi+h*(k1+2*k2+2*k3+k4)/6; outfileyi=yin; i+; while(i(b-a)/h+1);/此处为4阶Adams显式方法的通用程序void AB4(double a,double b,double h,double c) double x100,y100,y1100; dou

14、ble k1,k2,k3,k4; y0=c; x0=a; outfile输出4阶Adams显式方法的程序结果:n; for(int i=0;i=2;i+) xi=x0+i*h; k1=f1(xi,yi); k2=f1(xi+h/2),(yi+h*k1/2); k3=f1(xi+h/2),(yi+h*k2/2); k4=f1(xi+h),(yi+h*k3); yi+1=yi+h*(k1+2*k2+2*k3+k4)/6; int j=3; y10=y0; y11=y1; y12=y2; y13=y3; do xj=x0+j*h; y1j+1=y1j+(55*f1(xj,y1j)-59*f1(xj-

15、1,y1j-1)+37*f1(xj-2,y1j-2)-9*f1(xj-3,y1j-3)*h/24; outfiley1j=y1jn; j+; while(j(b-a)/h+1);/主函数void main(void) double a,b,h,c; coutabhc; accurate(a,b,h); RK4(a,b,h,c); AB4(a,b,h,c);float si(int u,int v) float sum=0; int q; for(q=0;qk;q+) sum+=auq*aqv; sum=auv-sum; return sum;void exchange(int g) int t

16、; float temp; for(t=0;tn;t+) temp=akt; akt=agt; agt=temp; void analyze() int t; float si(int u,int v); for(t=k;tn;t+) akt=si(k,t); for(t=(k+1);t-1;t-) sum=0; for(z=(t+1);zm;z+) sum+=atz*xz; xt=(float)(atm-sum)/att; （5）由经典Runge-Kutta公式得出的结果列在下面的表格中，以及精确值y(xi)和精确值和数值解的误差：ixiyiy(xi)|y(xi)-yi|0033010.12

17、.9972.9971.87138e-00720.22.976192.976193.91665e-00730.32.921132.921137.58342e-00740.42.819552.819551.61101e-00650.52.666662.66667 3.17735e-00660.62.46712.467115.00551e-00670.7 2.2338 2.23385.77233e-00680.81.984121.984134.12954e-00690.91.735111.735111.15554e-007101.01.500011.55.80668e-006111.11.28701

18、1.2871.13075e-005121.21.099721.099711.54242e-005131.30.9383970.938381.77272e-005141.40.80130.8012821.83754e-005151.50.6857320.6857141.78e-005由AB4方法得出的结果为：Y10=3 y11=2.997 y12=2.97619 y13=2.92113 y14=2.81839y15=2.66467 y16=2.4652 y17=2.23308 y18=1.98495 y19=1.73704y110=1.50219 y111=1.28876 y112=1.1007

19、2 y113=0.93871 y114=0.801135y115=0.685335（6）本次上机作业让我知道了在遇到复杂问题中，无法给出解析式的情况下，要学会灵活使用各种微分数值解法，并且能计算出不同方法的精度大小。第七章1）编制用Crank-Nicolson格式求抛物线方程2u/2x = f(x,t) (0x1, 0u(x,0) = (0u(0,t) = ，u（1,t）= (0t数值解的通用程序。2）就a=1, f(x,t)=0, =exp(x), =exp(t), =exp(1+t),M=40,N=40,输入点（0.2，1.0），（0.4，1.0），（0.6，1.0），（.8，1.0）4点

20、处u（x,t）的近似值。3） 已知所给方程的精确解为u(x,t)=exp(x+t),将步长反复二分，从（0.2，1.0），（0.4，1.0），（0.6，1.0），（0.8，1.0）4点处精确解与数值解的误差观察当步 长缩小一半时，误差以什么规律缩小。解：（1）#include #include float h=0.025,k=0.025;int m=40;int n=40; float y4040,r=a*k/(h*h);void Input() int i,j; coutLoading Input Data.endl; for(i=0;im;i+) for(j=0;jn;j+) if (i=

21、j) aij=1+r; for(j=0;jn;j+) if (j=i+1)|(i=j+1) aij=-r/2; int main() Input(); /read data int r,i; for(k=0;k(m-1);k+) int select(); /select main element r=select(); void exchange(int g); exchange(r); /exchange void analyze(); analyze(); /analyze void ret(); ret(); / replace back coutThe solution vector

22、 is below:endl; for(i=0;im;i+) coutxi=xiendl; return 0;int select() int f,t; float max; f=k; float si(int u,int v); max=float(fabs(si(k,k); for(t=(k+1);t(m-1);t+) if(maxfabs(si(t,k) max=float(fabs(si(t,k); f=t; return f;float si(int u,int v) float sum=0; int q; for(q=0;qk;q+) sum+=auq*aqv; sum=auv-s

23、um; return sum;void exchange(int g) int t; float temp; for(t=0;tn;t+) temp=akt; akt=agt; agt=temp; void analyze() int t; float si(int u,int v); for(t=k;tn;t+) akt=si(k,t); for(t=(k+1);t-1;t-) sum=0; for(z=(t+1);zm;z+) sum+=atz*xz; xt=(float)(atm-sum)/att; （2）结果：u(x,y)在所要求4点上的近似值：3.3201482664.0552499874.9530861226.049686221（3）精确解与数值解的误差随步长减小的变化情况：0.00050076540.00079890400.00085769590.00061934780.00012533150.00020001270.00021471660.00015497690.00003134340.00005002030.00005369740.00003875700.00000783650.00001250610.00001342