东南大学数值分析上机实验题(下).doc

《东南大学数值分析上机实验题(下).doc》由会员分享，可在线阅读，更多相关《东南大学数值分析上机实验题(下).doc（13页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

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 c99,e99,h99,u99

2、,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-x1j-1)*(x-x1j-1)*(x

3、-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); for (i=1;in;i+)

4、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;in;i+) aii-1=ui; a

5、ii+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; coutS(x)=; printY(f1i)

6、; 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+) if (i+0.5=x10) ca

7、lculation(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.0404019(x-2)(x-2)-0

8、.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.149051(x-6)-0.657537(x

9、-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 S(1.5)=3.67843 S (2.

10、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）中的四种方法进行计算，并将计算结果和精确解作比较；（6）通过本上机

11、题，你能得到哪些结论？解：#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(double a,double b,double h)d

12、ouble x100,accurate100;x0=a;int i=0;outfile输出函数准确值的程序结果:n;doxi=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;outfile输出经典Runge-Kutta公式的程序结果:n;doxi=x0+i

13、*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;double k1,k2,k3,k4;y0=c;x0=a;outfile输出4阶Adams显式方法的程序结果:n;

14、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;doxj=x0+j*h;y1j+1=y1j+(55*f1(xj,y1j)-59*f1(xj-1,y1j-1)+37*f1(xj-2,y1j-2)-9*f1(xj-3,y1j-3)*h/24;outfiley1j=y1jn;j+;whil

15、e(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; float temp;for(t=0;tn;t+)temp=akt;akt=agt;agt=temp;void analyze()int t;float si(int

16、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.9972.9971.87138e-00720.22.976192.976193.91665e-00730.32.921132.921137.58342e-00740.42.819552.819551

17、.61101e-00650.52.666662.666673.17735e-00660.62.46712.467115.00551e-00670.7 2.23382.23385.77233e-00680.81.984121.984134.12954e-00690.91.735111.735111.15554e-007101.01.500011.55.80668e-006111.11.287011.2871.13075e-005121.21.099721.099711.54242e-005131.30.9383970.938381.77272e-005141.40.80130.8012821.8

18、3754e-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.10072 y113=0.93871 y114=0.801135y115=0.685335（6） 本次上机作业让我知道了在遇到复杂问题中，无法给出解析式的情况下，要学会灵活使用各种微分数值解法，并且能计算出不同方

19、法的精度大小。第七章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点处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点处精确解与数值解的误差观

20、察当步 长缩小一半时，误差以什么规律缩小。解：（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=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

21、(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 is below:endl; for(i=0;im;i+) coutxi=xiendl; return 0;int select() int f,t; float max; f=k; float si(int

22、 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-sum; return sum;void exchange(int g) int t; float temp; for(t=0;tn;t+) temp=akt; akt=agt; agt=temp; void a

23、nalyze() 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.00001342550.00000969010.00000195920.00000312660.00000335650.00000242260.00000048980.00000078170.00000083910.00000060560.00000012240.00000019540.00000020980.0000001514