东南大学数值分析上机实验题下.docx
《东南大学数值分析上机实验题下.docx》由会员分享，可在线阅读，更多相关《东南大学数值分析上机实验题下.docx（17页珍藏版）》请在冰豆网上搜索。
东南大学数值分析上机实验题下
数值分析上机报告
姓名：
学号：
2013年12月22日
第四章
38.（上机题）3次样条插值函数
（1）编制求第一型3次样条插值函数的通用程序;
（2）已知汽车曲线型值点的数据如下：
0
1
2
3
4
5
6
7
8
9
10
2.51
3.30
4.04
4.70
5.22
5.54
5.78
5.40
5.57
5.70
5.80
端点条件为
=0.8，
=0.2。
用所编制程序求车门的3次样条插值函数S（x）,并打印出S（i+0.5）（i=0,1,…9）。
解：
（1）
#include
#include
floatx1[100],f1[100],f2[99],f3[98],m[100],a[100][101],x,d[100];
floatc[99],e[99],h[99],u[99],w[99],y_0,y_n,arr,s;
inti,j,k,n,q;
voidselectprint（floaty）
{
if（（y>0）&&（y!
=1））cout<<"+"<elseif（y==1）cout<<"+";
elseif（y<0）cout<}
voidprintY（floaty）{
if（y!
=0）cout<}
floatcalculation（floatx）{
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;
}
returns;
}
voidmain（）
{
do{
cout<<"请输入n值:
";
cin>>n;
if（（n>100）||（n<1））cout<<"请重新输入整数（1..100）:
"<}
while（（n>100）||（n<1））;
cout<<"请输入Xi（i=0,1,...,"<";
for（i=0;i<=n;i++）cin>>x1[i];
cout<cout<<"请输入Yi（i=0,1,...,"<";
for（i=0;i<=n;i++）cin>>f1[i];
cout<cout<<"请输入第一型边界条件Y'0，Y'n:
";
cin>>y_0>>y_n;
cout<for（i=0;ifor（i=1;ifor（i=1;ifor（i=0;ifor（i=0;ifor（i=1;id[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;
elsea[i][j]=0;
a[0][1]=1;
a[n][n-1]=1;
for（i=1;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;ic[i]=（float）（f2[i]-（1.0/3.0*m[i]+1.0/6.0*m[i+1]）*h[i]）;
for（i=0;ie[i]=（m[i+1]-m[i]）/（6*h[i]）;
for（i=0;i{
cout<<"X属于区间["<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<}
cout<<"针对本题，计算S（i+0.5）（i=0,1,...,9）:
"<for（i=0;i<10;i++）
{
if（（i+0.5<=x1[n]）&&（i+0.5>=x1[0]））
{
calculation（（float）（i+0.5））;
cout<<"S（"<}
elsecout<};
cout<}
（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.90856S（1.5）=3.67843S（2.5）=4.38147
S（3.5）=4.98819S（4.5）=5.38328S（5.5）=5.7237
S（6.5）=5.59441S（7.5）=5.42989S（8.5）=5.65976
S（9.5）=5.7323
第六章
21．（上机题）常微分方程初值问题数值解
（1）编制RK4方法的通用程序；
（2）编制AB4方法的通用程序（由RK4提供初值）；
（3）编制AB4-AM4预测校正方法的通用程序（由RK4提供初值）；
（4）编制带改进的AB4-AM4预测校正方法的通用程序（由RK4提供初值）;
（5）对于初值问题
取步长
应用
（1）~（4）中的四种方法进行计算，并将计算结果和精确解
作比较；
（6）通过本上机题，你能得到哪些结论？
解：
#include
#include
#include
#include
ofstreamoutfile（"data.txt"）;
//此处定义函数f（x,y）的表达式
//用户可以自己设定所需要求得函数表达式
doublef1（doublex,doubley）
{
doublef1;
f1=（-1）*x*x*y*y;
returnf1;
}
//此处定义求函数精确解的函数表达式
doublef2（doublex）
{
doublef2;
f2=3/（1+x*x*x）;
returnf2;
}
//此处为精确求函数解的通用程序
voidaccurate（doublea,doubleb,doubleh）
{
doublex[100],accurate[100];
x[0]=a;
inti=0;
outfile<<"输出函数准确值的程序结果:
\n";
do{
x[i]=x[0]+i*h;
accurate[i]=f2（x[i]）;
outfile<<"accurate["<i++;
}while（i<（b-a）/h+1）;
}
//此处为经典Runge-Kutta公式的通用程序
voidRK4（doublea,doubleb,doubleh,doublec）
{
inti=0;
doublek1,k2,k3,k4;
doublex[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++;
}while（i<（b-a）/h+1）;
}
//此处为4阶Adams显式方法的通用程序
voidAB4（doublea,doubleb,doubleh,doublec）
{
doublex[100],y[100],y1[100];
doublek1,k2,k3,k4;
y[0]=c;
x[0]=a;
outfile<<"输出4阶Adams显式方法的程序结果:
\n";
for（inti=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;
}
intj=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++;
}while（j<（b-a）/h+1）;
}
//主函数
voidmain（void）
{
doublea,b,h,c;
cout<<"输入上下区间、步长和初始值:
\n";
cin>>a>>b>>h>>c;
accurate（a,b,h）;
RK4（a,b,h,c）;
AB4（a,b,h,c）;
}
floatsi（intu,intv）
{
floatsum=0;intq;
for（q=0;qsum+=a[u][q]*a[q][v];
sum=a[u][v]-sum;
returnsum;
}
voidexchange（intg）
{
intt;floattemp;
for（t=0;t{
temp=a[k][t];
a[k][t]=a[g][t];
a[g][t]=temp;
}
}
voidanalyze（）
{
intt;
floatsi（intu,intv）;
for（t=k;ta[k][t]=si（k,t）;
for（t=（k+1）;ta[t][k]=（float）（si（t,k）/a[k][k]）;
}
voidret（）
{
intt,z;floatsum;
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）;zsum+=a[t][z]*x[z];
x[t]=（float）（a[t][m]-sum）/a[t][t];
}
}
（5）
由经典Runge-Kutta公式得出的结果列在下面的表格中，以及精确值y（xi）和精确值和数值解的误差：
i
xi
yi
y（xi）
|y（xi）-yi|
0
0
3
3
0
1
0.1
2.997
2.997
1.87138e-007
2
0.2
2.97619
2.97619
3.91665e-007
3
0.3
2.92113
2.92113
7.58342e-007
4
0.4
2.81955
2.81955
1.61101e-006
5
0.5
2.66666
2.66667
3.17735e-006
6
0.6
2.4671
2.46711
5.00551e-006
7
0.7
2.2338
2.2338
5.77233e-006
8
0.8
1.98412
1.98413
4.12954e-006
9
0.9
1.73511
1.73511
1.15554e-007
10
1.0
1.50001
1.5
5.80668e-006
11
1.1
1.28701
1.287
1.13075e-005
12
1.2
1.09972
1.09971
1.54242e-005
13
1.3
0.938397
0.93838
1.77272e-005
14
1.4
0.8013
0.801282
1.83754e-005
15
1.5
0.685732
0.685714
1.78e-005
由AB4方法得出的结果为：
Y1[0]=3y1[1]=2.997y1[2]=2.97619y1[3]=2.92113y1[4]=2.81839
y1[5]=2.66467y1[6]=2.4652y1[7]=2.23308y1[8]=1.98495y1[9]=1.73704
y1[10]=1.50219y1[11]=1.28876y1[12]=1.10072y1[13]=0.93871y1[14]=0.801135
y1[15]=0.685335
（6）本次上机作业让我知道了在遇到复杂问题中，无法给出解析式的情况下，要学会灵活使用各种微分数值解法，并且能计算出不同方法的精度大小。
第七章
1）编制用Crank-Nicolson格式求抛物线方程
2u/
2x=f（x,t）（0u（x,0）=
（0
u（0,t）=
，u（1,t）=
（0数值解的通用程序。
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点处精确解与数值解的误差观察当步长缩小一半时，误差以什么规律缩小。
解：
（1）#include
#include
floath=0.025,k=0.025;
intm=40;intn=40;
floaty[40][40],r=a*k/（h*h）;
voidInput（）
{
inti,j;
cout<<"LoadingInputData..."<for（i=0;i{
for（j=0;jif（i==j）a[i][j]=1+r;
for（j=0;jif（（j=i+1）||（i=j+1））a[i][j]=-r/2;
}
}
intmain（）
{
Input（）;//readdata
intr,i;
for（k=0;k<（m-1）;k++）
{
intselect（）;//selectmainelement
r=select（）;
voidexchange（intg）;
exchange（r）;//exchange
voidanalyze（）;
analyze（）;//analyze
}
voidret（）;
ret（）;//replaceback
cout<<"Thesolutionvectorisbelow:
"<for（i=0;icout<<"x["<return0;
}
intselect（）
{
intf,t;floatmax;
f=k;
floatsi（intu,intv）;
max=float（fabs（si（k,k）））;
for（t=（k+1）;t<（m-1）;t++）
{
if（max{
max=float（fabs（si（t,k）））;
f=t;
}
}
returnf;
}
floatsi（intu,intv）
{
floatsum=0;intq;
for（q=0;qsum+=a[u][q]*a[q][v];
sum=a[u][v]-sum;
returnsum;
}
voidexchange（intg）
{
intt;floattemp;
for（t=0;t{
temp=a[k][t];
a[k][t]=a[g][t];
a[g][t]=temp;
}
}
voidanalyze（）
{
intt;
floatsi（intu,intv）;
for（t=k;ta[k][t]=si（k,t）;
for（t=（k+1）;ta[t][k]=（float）（si（t,k）/a[k][k]）;
}
voidret（）
{
intt,z;floatsum;
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）;zsum+=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.00001342
升级会员 