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

数值分析上机报告

姓名： 学号：

2013年12月22日

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

第四章

38.（上机题）3次样条插值函数

（1）编制求第一型3次样条插值函数的通用程序;

''

端点条件为y0=0.8，y10=0.2。用所编制程序求车门的3次样条插值函数S(x),并打印出

S(i+0.5)(i=0,1,…9)。

解：

（1）

#include <iostream.h> #include <math.h>

float x1[100],f1[100],f2[99],f3[98],m[100],a[100][101],x,d[100]; float c[99],e[99],h[99],u[99],w[99],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;j<=n;j++) if (x<=x1[j]) {

s=(float)(f1[j-1]+c[j-1]*(x-x1[j-1])+m[j-1]/2.0*(x-x1[j-1])*(x-x1[j-1])+e[j-1]*(x-x1[j-1])*(x-x1[j-1])*(x-x1[j-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));

cout<<"请输入Xi(i=0,1,...,"<<n<<"):"; for (i=0;i<=n;i++) cin>>x1[i]; cout<<endl;

cout<<"请输入Yi(i=0,1,...,"<<n<<"n):"; for (i=0;i<=n;i++) cin>>f1[i]; cout<<endl;

cout<<"请输入第一型边界条件Y'0，Y'n:"; cin>>y_0>>y_n; cout<<endl;

for (i=0;i<n;i++) h[i]=x1[i+1]-x1[i];

for (i=1;i<n;i++) u[i]=(float) (h[i-1]/(h[i-1]+h[i])); for (i=1;i<n;i++) w[i]=(float) (1.0-u[i]);

for (i=0;i<n;i++) f2[i]=(f1[i+1]-f1[i])/h[i]; //一阶差商 for (i=0;i<n-1;i++) f3[i]=(f2[i+1]-f2[i])/(x1[i+2]-x1[i]); //二阶差商 for (i=1;i<n;i++) d[i]=6*f3[i-1]; //求出所有的d[i] d[0]=6*(f2[0]-y_0)/h[0]; d[n]=6*(y_n-f2[n-1])/h[n-1]; for (i=0;i<=n;i++) for (j=0;j<=n;j++)

if (i==j) a[i][j]=2; else a[i][j]=0; a[0][1]=1; a[n][n-1]=1;

for (i=1;i<n;i++) {

a[i][i-1]=u[i]; a[i][i+1]=w[i]; }

for (i=0;i<=n;i++) a[i][n+1]=d[i];

for (i=1;i<=n;i++) //用追赶法解方程,得M[i] {

arr=a[i][i-1];

for (j=0;j<=n+1;j++)

a[i][j]=a[i][j]-a[i-1][j]*arr/a[i-1][i-1]; }

m[n]=a[n][n+1]/a[n][n];

for (i=n-1;i>=0;i--) m[i]=(a[i][n+1]-a[i][i+1]*m[i+1])/a[i][i];

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

for (i=0;i<n;i++) //计算S（x）的表达式

c[i]=(float) (f2[i]-(1.0/3.0*m[i]+1.0/6.0*m[i+1])*h[i]); for (i=0;i<n;i++)

e[i]=(m[i+1]-m[i])/(6*h[i]); for (i=0;i<n;i++) {

cout<<"X属于区间["<<x1[i]<<","<<x1[i+1]<<"]时"<<endl<<endl; cout<<"S(x)="; printY(f1[i]); if (c[i]!=0) {

selectprint(c[i]); cout<<"(x"; printY(-x1[i]); cout<<")"; }

if (m[i]!=0) {

selectprint((float)(m[i]/2.0)); for (q=0;q<2;q++) {

cout<<"(x"; printY(-x1[i]); cout<<")"; } }

if (e[i]!=0) {

selectprint(e[i]); for (q=0;q<3;q++) {

cout<<"(x"; printY(-x1[i]); cout<<")"; } }

cout<<endl<<endl; }

cout<<"针对本题，计算S(i+0.5)(i=0,1,...,9):"<<endl; for (i=0;i<10;i++) {

if ((i+0.5<=x1[n])&&(i+0.5>=x1[0])) {

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

calculation((float)(i+0.5));

cout<<"S("<<i+0.5<<")="<<s<<endl; }

else cout<<i+0.5<<"超出定义域"<<endl; }; cout<<endl; } （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.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-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.5)=4.38147 S(3.5)=4.98819 S(4.5)=5.38328 S(5.5)=5.7237 S(6.5)=5.59441 S(7.5)=5.42989 S(8.5)=5.65976 S(9.5)=5.7323

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

第六章

21．（上机题）常微分方程初值问题数值解 （1）编制RK4方法的通用程序；

（2）编制AB4方法的通用程序（由RK4提供初值）；

（3）编制AB4-AM4预测校正方法的通用程序（由RK4提供初值）；

（4）编制带改进的AB4-AM4预测校正方法的通用程序（由RK4提供初值）; （5）对于初值问题

取步长h 0.1,应用（1）~（4）中的四种方法进行计算，并将计算结果和精确解

y(x) 3/(1 x3)作比较；

（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=(-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 x[100],accurate[100]; x[0]=a; int i=0; outfile<<"输出函数准确值的程序结果:\n"; do{ x[i]=x[0]+i*h; accurate[i]=f2(x[i]); outfile<<"accurate["<<i<<"]="<<accurate[i]<<'\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 x[100],y[100]; y[0]=c; x[0]=a; outfile<<"输出经典Runge-Kutta公式的程序结果:\n"; do { x[i]=x[0]+i*h; k1=f1(x[i],y[i]); k2=f1((x[i]+h/2),(y[i]+h*k1/2)); k3=f1((x[i]+h/2),(y[i]+h*k2/2)); k4=f1((x[i]+h),(y[i]+h*k3)); y[i+1]=y[i]+h*(k1+2*k2+2*k3+k4)/6; outfile<<"y"<<"["<<i<<"]="<<y[i]<<'\n'; i++; }while(i<(b-a)/h+1); }

//此处为4阶Adams显式方法的通用程序 void AB4(double a,double b,double h,double c) { double x[100],y[100],y1[100]; double k1,k2,k3,k4; y[0]=c; x[0]=a; outfile<<"输出4阶Adams显式方法的程序结果:\n"; for(int i=0;i<=2;i++)

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

{ x[i]=x[0]+i*h; k1=f1(x[i],y[i]); k2=f1((x[i]+h/2),(y[i]+h*k1/2)); k3=f1((x[i]+h/2),(y[i]+h*k2/2)); k4=f1((x[i]+h),(y[i]+h*k3)); y[i+1]=y[i]+h*(k1+2*k2+2*k3+k4)/6; } int j=3; y1[0]=y[0]; y1[1]=y[1]; y1[2]=y[2]; y1[3]=y[3]; do { x[j]=x[0]+j*h; y1[j+1]=y1[j]+(55*f1(x[j],y1[j])-59*f1(x[j-1],y1[j-1])+37*f1(x[j-2],y1[j-2])-9*f1(x[j-3],y1[j-3]))*h/24; outfile<<"y1"<<"["<<j<<"]="<<y1[j]<<'\n'; j++; }while(j<(b-a)/h+1); }

//主函数

void main(void) { double a,b,h,c; cout<<"输入上下区间、步长和初始值:\n"; cin>>a>>b>>h>>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+=a[u][q]*a[q][v]; sum=a[u][v]-sum; return sum; }

void exchange(int g) {

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

int t; float temp; for(t=0;t<n;t++) { temp=a[k][t]; a[k][t]=a[g][t]; a[g][t]=temp; } }

void analyze() { int t; float si(int u,int v); for(t=k;t<n;t++) a[k][t]=si(k,t); for(t=(k+1);t<m;t++) a[t][k]=(float)(si(t,k)/a[k][k]); }

void ret() { int t,z;float sum; x[m-1]=(float)a[m-1][m]/a[m-1][m-1]; for(t=(m-2);t>-1;t--) { sum=0; for(z=(t+1);z<m;z++) sum+=a[t][z]*x[z]; x[t]=(float)(a[t][m]-sum)/a[t][t]; } }

（5）

由经典Runge-Kutta公式得出的结果列在下面的表格中，以及精确值y(xi)和精确值和数值解的误差：

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

由AB4方法得出的结果为：

Y1[0]=3 y1[1]=2.997 y1[2]=2.97619 y1[3]=2.92113 y1[4]=2.81839 y1[5]=2.66467 y1[6]=2.4652 y1[7]=2.23308 y1[8]=1.98495 y1[9]=1.73704 y1[10]=1.50219 y1[11]=1.28876 y1[12]=1.10072 y1[13]=0.93871 y1[14]=0.801135 y1[15]=0.685335

（6）本次上机作业让我知道了在遇到复杂问题中，无法给出解析式的情况下，要学会灵活使用各种微分数值解法，并且能计算出不同方法的精度大小。

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

第七章

1）编制用Crank-Nicolson格式求抛物线方程

u/ t - a 2u/ 2x = f(x,t) (0<x<1, 0 t T)

u(x,0) = (x) (0 x 1) u(0,t) = x ，u（1,t）= t (0<t 1) 数值解的通用程序。

2）就a=1, f(x,t)=0, x =exp(x), t =exp(t), 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点处精确解与数值解的误差观察当步 长缩小一半时，误差以什么规律缩小。 解：

（1）#include <iostream> #include <math.h>

float h=0.025,k=0.025; int m=40;int n=40;

float y[40][40],r=a*k/(h*h);

void Input() {

int i,j;

cout<<"Loading Input Data..."<<endl; for(i=0;i<m;i++) {

for(j=0;j<n;j++)

if (i==j) a[i][j]=1+r; for(j=0;j<n;j++)

if ((j=i+1)||(i=j+1)) a[i][j]=-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

cout<<"The solution vector is below:"<<endl; for(i=0;i<m;i++)

cout<<"x["<<i<<"]="<<x[i]<<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))); f=t; } } return f; }

float si(int u,int v) {

float sum=0; int q; for(q=0;q<k;q++)

sum+=a[u][q]*a[q][v]; sum=a[u][v]-sum; return sum; }

void exchange(int g)

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

{

int t; float temp; for(t=0;t<n;t++) {

temp=a[k][t]; a[k][t]=a[g][t]; a[g][t]=temp; } }

void analyze() {

int t;

float si(int u,int v); for(t=k;t<n;t++) a[k][t]=si(k,t); for(t=(k+1);t<m;t++)

a[t][k]=(float)(si(t,k)/a[k][k]); }

void ret() {

int t,z;float sum;

x[m-1]=(float)a[m-1][m]/a[m-1][m-1]; for(t=(m-2);t>-1;t--) {

sum=0;

for(z=(t+1);z<m;z++) sum+=a[t][z]*x[z];

x[t]=(float)(a[t][m]-sum)/a[t][t]; } }

（2）结果：

u(x,y)在所要求4点上的近似值：

3.320148266 4.055249987 4.953086122 6.049686221 （3）精确解与数值解的误差随步长减小的变化情况： 0.0005007654 0.0007989040 0.0008576959 0.0006193478 0.0001253315 0.0002000127 0.0002147166 0.0001549769 0.0000313434 0.0000500203 0.0000536974 0.0000387570 0.0000078365 0.0000125061 0.0000134255 0.0000096901 0.0000019592 0.0000031266 0.0000033565 0.0000024226 0.0000004898 0.0000007817 0.0000008391 0.0000006056 0.0000001224 0.0000001954 0.0000002098 0.0000001514