東南大學數值分析上機實驗題下

1、. 數值分析上機報告XX：學號： 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 <iostream.h>*include <math.h>floatx1100,f1100,f299,f398,m100,a100101

2、,x,d100;float c99,e99,h99,u99,w99,y_0,y_n,arr,s;int i,j,k,n,q;void selectprint(float y)if (y>0)&&(y!=1) cout<<"+"<<y; else if (y=1) cout<<"+" else if (y<0) cout<<y;void printY(float y)if (y!=0) cout<<y;float calculation(float x)for (j=1

3、;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-x1j-1); break; return s;void main()do cout<<"請輸入n值:" cin>>n; if (n>100)|(n<1) cout<<"請重新輸入整數(1.100):"<<endl; while (n>100)|(n<1)

4、; cout<<"請輸入Xi(i=0,1,.,"<<n<<"):" for (i=0;i<=n;i+) cin>>x1i; cout<<endl; cout<<"請輸入Yi(i=0,1,.,"<<n<<"n):" for (i=0;i<=n;i+) cin>>f1i; cout<<endl; cout<<"請輸入第一型邊界條件Y'0，Y'n:&qu

5、ot; cin>>y_0>>y_n; cout<<endl; for (i=0;i<n;i+) hi=x1i+1-x1i; for (i=1;i<n;i+) ui=(float) (hi-1/(hi-1+hi); for (i=1;i<n;i+) wi=(float) (1.0-ui); for (i=0;i<n;i+) f2i=(f1i+1-f1i)/hi; /一階差商 for (i=0;i<n-1;i+) f3i=(f2i+1-f2i)/(x1i+2-x1i); /二階差商 for (i=1;i<n;i+) di=6*

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;i<n;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<=n+1;j+) aij=aij-ai-1j*arr

7、/ai-1i-1; mn=ann+1/ann; for (i=n-1;i>=0;i-) mi=(ain+1-aii+1*mi+1)/aii; for (i=0;i<n;i+) /計算S（x）的表達式 ci=(float) (f2i-(1.0/3.0*mi+1.0/6.0*mi+1)*hi); for (i=0;i<n;i+) ei=(mi+1-mi)/(6*hi); for (i=0;i<n;i+) cout<<"X屬于區(qū)間"<<x1i<<","<<x1i+1<<&quo

8、t;時"<<endl<<endl; cout<<"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;q<2;q+) cout<<"(x" printY(-x1i); cout<<")" i

9、f (ei!=0) selectprint(ei); for (q=0;q<3;q+) cout<<"(x" printY(-x1i); cout<<")" cout<<endl<<endl; cout<<"針對本題，計算S(i+0.5)(i=0,1,.,9):"<<endl;for (i=0;i<10;i+) if (i+0.5<=x1n)&&(i+0.5>=x10) calculation(float)(i+0.5);

10、cout<<"S("<<i+0.5<<")="<<s<<endl; else cout<<i+0.5<<"超出定義域"<<endl;cout<<endl;（2）編制的程序求車門的3次樣條插值函數S(x)：x屬于區(qū)間0,1時；S(x)=2.51+0.8(x)-0.0014861(x)(x)-0.00851395(x)(x)(x)x屬于區(qū)間1,2時；S(x)=3.3+0.771486(x-1)-0.027028(x-1)(x-1)-

11、0.00445799(x-1)(x-1)(x-1)x屬于區(qū)間2,3時；S(x)=4.04+0.704056(x-2)-0.0404019(x-2)(x-2)-0.0036543(x-2)(x-2)(x-2)x屬于區(qū)間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屬于區(qū)間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屬于區(qū)間5,6時；S(x)=5.54+0.360567(x-5)+0.147919(x

12、-5)(x-5)-0.268485(x-5)(x-5)(x-5)x屬于區(qū)間6,7時；S(x)=5.78-0.149051(x-6)-0.657537(x-6)(x-6)+0.426588(x-6)(x-6)(x-6)x屬于區(qū)間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屬于區(qū)間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屬于區(qū)間9,10時；S(x)=5.7+0.058376(x-9)-0.0167

13、508(x-9)(x-9)+0.0583752(x-9)(x-9)(x-9)S(0.5)=2.90856 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預測校

14、正方法的通用程序（由RK4提供初值）;（5）對于初值問題取步長,應用（1）（4）中的四種方法進行計算，并將計算結果和精確解作比較；（6）通過本上機題，你能得到哪些結論.解：*include<iostream.h>*include<fstream.h>*include<stdlib.h>*include<math.h>ofstream outfile("data.txt");/此處定義函數f(x,y)的表達式/用戶可以自己設定所需要求得函數表達式double f1(double x,double y)double f1;f1=(

15、-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)double x100,accurate100;x0=a;int i=0;outfile<<"輸出函數準確值的程序結果:n"doxi=x0+i*h;accuratei=f2(xi);outfile<<"accurate"<<i

16、<<"="<<accuratei<<'n'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*h;k1=f1(xi,yi);k2=f1(xi+h/2),(yi+h

17、*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;outfile<<"y"<<""<<i<<"="<<yi<<'n'i+;while(i<(b-a)/h+1);/此處為4階Adams顯式方法的通用程序void AB4(double a,double b,double h,double c)double x100,y100,y

18、1100;double 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;doxj=x0+j*h;y1j+1=y1j+(55*f1(xj,y

19、1j)-59*f1(xj-1,y1j-1)+37*f1(xj-2,y1j-2)-9*f1(xj-3,y1j-3)*h/24;outfile<<"y1"<<""<<j<<"="<<y1j<<'n'j+;while(j<(b-a)/h+1);/主函數void main(void)double a,b,h,c;cout<<"輸入上下區(qū)間、步長和初始值:n"cin>>a>>b>>h&

20、gt;>c;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;q<k;q+)sum+=auq*aqv;sum=auv-sum;return sum;void exchange(int g)int t; float temp;for(t=0;t<n;t+)temp=akt;akt=agt;agt=temp;void analyze()int t;float si(int u,int v);for(t=k;t<n;t+)akt=si(k,t);f

21、or(t=(k+1);t<m;t+)atk=(float)(si(t,k)/akk);void ret()int t,z;float sum;xm-1=(float)am-1m/am-1m-1;for(t=(m-2);t>-1;t-)sum=0;for(z=(t+1);z<m;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.

22、976192.976193.91665e-00730.32.921132.921137.58342e-00740.42.819552.819551.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.09972

23、1.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.10072 y113=0.93871 y114=0.801135y

24、115=0.685335（6） 本次上機作業(yè)讓我知道了在遇到復雜問題中，無法給出解析式的情況下，要學會靈活使用各種微分數值解法，并且能計算出不同方法的精度大小。第七章1）編制用Crank-Nicolson格式求拋物線方程2u/2x = f(x,t) (0<x<1, 0u(x,0) = (0u(0,t) = ，u（1,t）= (0<t數值解的通用程序。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） 已知所給

25、方程的精確解為u(x,t)=exp(x+t),將步長反復二分，從（0.2，1.0），（0.4，1.0），（0.6，1.0），（0.8，1.0）4點處精確解與數值解的誤差觀察當步 長縮小一半時，誤差以什么規(guī)律縮小。解：（1）*include <iostream>*include <math.h>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; cout<<"Loading Input Data."<<end

26、l; for(i=0;i<m;i+) for(j=0;j<n;j+) if (i=j) aij=1+r; for(j=0;j<n;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 voi

27、d ret(); ret(); / replace back cout<<"The solution vector is below:"<<endl; for(i=0;i<m;i+) cout<<"x"<<i<<"="<<xi<<endl; 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(max<fabs(si(t,k) max=float(fabs(si(t,k);