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

数值分析上机练习

（以VC++6.0为操作平台）

第四章（4.38） 程序如下：

#include<iostream.h> void main(void) {

float x[11];//存放数组x[j] float y[11];//存放数组y[j] float h[11];//存放数组h[j] float u[11];//存放数组u[j] float v[11];//存放数组v[j] float d[11];//存放数组d[j] float M[11];//存放数组M[j] float b[11];// 存放数组b[j]

float t[11],l[11],yy[11],s[4],aa1,aa2,aa3,aa4; float s1[10]; int i,j,n; float xx;

cout<<"请输入n的值:\n"; cin>>n;

cout<<"输入数组x:\n"; for(i=0;i<=n;i++)cin>>x[i]; cout<<"输入数组y:\n"; for(i=0;i<=n;i++)cin>>y[i]; //输入端点值 float df[2];

cout<<"输入两个端点值:\n";

for(i=0;i<2;i++)cin>>df[i];//求出h[j]的值 for(j=0;j<=n-1;j++) {h[j]=x[j+1]-x[j];

cout<<'h'<<'['<<j<<']'<<'='<<h[j]<<'\t'; }

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

cout<<endl;

//求出u[j]和v[j]的初值 v[0]=1; u[n]=1;

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

{u[j]=h[j-1]/(h[j-1]+h[j]); v[j]=h[j]/(h[j-1]+h[j]); }

//求出d[j]的值 for(j=1;j<n;j++)

{d[j]=6*((y[j+1]-y[j])/h[j]-(y[j]-y[j-1])/h[j-1])/(h[j]+h[j-1]);} d[0]=6*((y[1]-y[0])/h[0]-df[0])/h[0];

d[n]=6*(df[1]-(y[n]-y[n-1])/h[n-1])/h[n-1]; for(j=1;j<=n;j++) {

cout<<'u'<<'['<<j<<']'<<'='<<u[j]<<'\t'; }

cout<<endl;

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

cout<<'v'<<'['<<j<<']'<<'='<<v[j]<<'\t'; }

cout<<endl;

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

cout<<'d'<<'['<<j<<']'<<'='<<d[j]<<'\t'; }

cout<<endl;

//利用书本上的追赶法求解方程组 for(i=0;i<=n;i++) {b[i]=2;} cout<<endl; t[0]=b[0]; yy[0]=d[0];

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

//消元过程

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

l[i]=u[i]/t[i-1]; t[i]=b[i]-l[i]*v[i-1]; yy[i]=d[i]-l[i]*yy[i-1]; }

//回代过程

M[n]=yy[n]/t[n]; for(i=n-1;i>=0;i--) {

M[i]=(yy[i]-v[i]*M[i+1])/t[i]; }

//将M[j]的值输出 for(i=0;i<=n;i++)

{cout<<'M'<<'['<<i<<']'<<'='<<M[i]<<endl;} //输出插值多项式的系数 for(j=0;j<n;j++) {

s[0]=y[j];

s[1]=(y[j+1]-y[j])/h[j]-(M[j]/3+M[j+1]/6)*h[j]; s[2]=M[j]/2;

s[3]=(M[j+1]-M[j])/(6*h[j]);

cout<<"当x的值在区间"<<'x'<<'['<<j<<']'<<"到"<<'x'<<'['<<(j+1)<<']'<<"时,输出插值多项式的系数:\n"; for(int k=0;k<4;k++) {

cout<<'s'<<'['<<k<<']'<<'='<<s[k]<<'\t'; }

cout<<endl; } }

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

程序结果：详见附图4.38jpg

编制的程序求车门的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

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

第六章（6.21） 程序如下：

#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';

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

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); }

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;

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

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); }

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

程序结果：

经典Runge-Kutta公式得出的结果列在下面的表格中，以及精确

由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.5021 y1[11]=1.28876 y1[12]=1.10072 y1[13]=0.93871 y1[14]=0.801135 y1[15]=0.685335

通过本上机题我明白了各种求微分方程的数值方法，经典Runge-Kutta公式，AB4方法以及AB4-AM4预测校正方法求解公式的精度是不同的。其中经典Runge-Kutta公式的精度，四阶Adams显式方法（AB4）均具有4阶精度。

第八章（8.12） 程序如下：

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

double CN_f(double x,double t) //f(x,y) {

return 0; }

double CN_phi(double x) //phi(x)

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

{

return exp(x); }

double CN_alpha(double t) //alpha(t) {

return exp(t); }

double CN_beta(double t) //beta(t) {

return exp(1+t); }

void ZG(double a[],double b[],double c[],double d[],int n1,double Ans[]) //追赶法解三角方程 {

double *beta,*y,*L; beta=new double[n1]; L=new double[n1]; y=new double[n1];

beta[0]=b[0]; y[0]=d[0];

for(int i=1;i<n1;i++) {

L[i]=a[i]/beta[i-1];

beta[i]=b[i]-L[i]*c[i-1]; y[i]=d[i]-L[i]*y[i-1]; }

Ans[n1-1]=y[n1-1]/beta[n1-1]; for(i=n1-2;i>=0;i--) {

Ans[i]=(y[i]-c[i]*Ans[i+1])/beta[i]; }

delete[] beta;

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

delete[] L; delete[] y; }

void CN(double A,double h,int M,int N_T,double tau,double Ans[]) //Crank-Nicolson格式 {

double r;

r=A*tau/(h*h);

double *a,*b,*c,*d,*uk,*uk1; a=new double[M-1]; b=new double[M-1]; c=new double[M-2]; d=new double[M-1]; uk=new double[M]; uk1=new double[M]; b[0]=1+r; c[0]=-r/2; a[M-2]=-r/2; b[M-2]=1+r;

for(int i=1;i<=M-3;i++) {

a[i]=-r/2; b[i]=1+r; c[i]=-r/2; }

for(i=1;i<=M-1;i++) {

uk[i]=CN_phi(i*h); }

for(i=1;i<=M-1;i++) {

Ans[i]=uk[i]; }

for(int k=0;k<=N_T-1;k++)

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

d[0]=(1-r)*uk[1]+r/2*uk[2]

+tau*CN_f(1*h,(k+0.5)*tau)

+r/2*(CN_alpha(k*tau)+CN_alpha((k+1)*tau)); d[M-2]=r/2*uk[M-2]+(1-r)*uk[M-1] +tau*CN_f((M-1)*h,(k+0.5)*tau)

+r/2*(CN_beta(k*tau)+CN_beta((k+1)*tau)); for(i=1;i<=M-3;i++) {

d[i]=r/2*uk[i]+(1-r)*uk[i+1]+r/2*uk[i+2] +tau*CN_f((i+1)*h,(k+0.5)*tau); }

ZG(a,b,c,d,M-1,uk1); for(i=1;i<=M-1;i++) {

Ans[(k+1)*(M+1)+i]=uk1[i-1]; }

for(i=1;i<=M-1;i++) {

uk[i]=uk1[i-1]; } }

for(k=0;k<=N_T;k++) {

Ans[k*(M+1)]=CN_alpha(k*tau); Ans[k*(M+1)+M]=CN_beta(k*tau); }

delete[] a; delete[] b; delete[] c; delete[] d; delete[] uk; delete[] uk1;

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

void main(void) {

double A,h,tau,*Ans; int M,N_T; A=1; M=40; N_T=40; h=1.0/M; tau=1.0/N_T;

Ans=new double[(M+1)*(N_T+1)]; CN(A,h,M,N_T,tau,Ans);

cout<<"u(x,y)在所要求4点上的近似值："<<endl; cout<<setprecision(10)<<Ans[40*(M+1)+8]<<setw(15); cout<<Ans[40*(M+1)+16]<<setw(15); cout<<Ans[40*(M+1)+24]<<setw(15);

cout<<Ans[40*(M+1)+32]<<setw(15)<<endl; delete[] Ans;

cout<<"精确解与数值解的误差随步长减小的变化情况："<<endl; N_T=5; M=5; int k=1;

for(int i=2;(M*=2)<1000;) {

N_T*=2; h=1.0/M; tau=1.0/N_T;

Ans=new double[(M+1)*(N_T+1)]; CN(A,h,M,N_T,tau,Ans); k*=2;

cout<<setiosflags(ios::fixed)

<<Ans[N_T*(M+1)+k]-exp(0.2+1.0)<<setw(15);

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

cout<<Ans[N_T*(M+1)+2*k]-exp(0.4+1.0)<<setw(15); cout<<Ans[N_T*(M+1)+3*k]-exp(0.6+1.0)<<setw(15); cout<<Ans[N_T*(M+1)+4*k]-exp(0.8+1.0)<<endl; delete[] Ans; } }

程序结果：

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

3.320148266 4.055249987 4.953086122 精确解与数值解的误差随步长减小的变化情况：0.0005007654 0.0007989040 0.0008576959 0.0001253315 0.0002000127 0.0002147166 0.0000313434 0.0000500203 0.0000536974 0.0000078365 0.0000125061 0.0000134255 0.0000019592 0.0000031266 0.0000033565 0.0000004898 0.0000007817 0.0000008391 0.0000001224 0.0000001954 0.0000002098 Press any key to continue

6.049686221 0.0006193478 0.0001549769 0.0000387570 0.0000096901 0.0000024226 0.0000006056 0.0000001514