13级数值分析实验报告.docx
- 文档编号:13773741
- 上传时间:2023-06-17
- 格式:DOCX
- 页数:29
- 大小:316.73KB
13级数值分析实验报告.docx
《13级数值分析实验报告.docx》由会员分享,可在线阅读,更多相关《13级数值分析实验报告.docx(29页珍藏版)》请在冰点文库上搜索。
13级数值分析实验报告
数值分析上机实验报告
课程名称:
学院:
专业:
姓名:
学号:
年级:
任课教师:
一,用Newton法求方程:
x7-28x4+14=0在(0.1,1.9)中的近似根(初始近似值取为区间端点,迭代6次或误差小于0.00001.)
解:
(1)解题方法及理论依据:
Newton迭代法定理:
设函数在有限区间
上二阶导数存在,且满足条件:
;
在区间
上不变号;
;
其中
是
中使
达到的一个。
则对任意初始近似值
,由Newton迭代过程
所生的迭代序列
平方收敛于方程
在区间
上的唯一解
。
(2)计算机程序(C语言):
#include
#include
voidmain()
{intflag=1,c;//设置flag控制while循环(即控制整个程序是否继续运行)
while(flag==1)
{floatx,y,f,f1;
inti=0;
printf("\n请输入x的值:
\n");
scanf("%f",&x);//输入x的初值
for(;(f1!
=0)&&(fabs(y-x)>=0.000001);)//迭代循环控制条件
{y=x;
f=pow(x,7)-28*pow(x,4)+14;//原方程
f1=7*pow(x,6)-28*4*pow(x,3);//原方程的导数方程
x=x-f/f1;//牛顿迭代公式
i++;//累计迭代的次数
printf("x=%f,y=%f\n",x,y);//输出每次迭代的值
}
printf("\n方程的近似根:
y=%f\n",y);
printf("\n迭代次数:
%d\n",i);
printf("是否继续(Y/N)?
");//提示是否输入下一个x的初值
getchar();
c=getchar();
if(c=='N'||c=='n')//如果输入N或n,则flag=0退出程序;则继续运行程序
flag=0;
elseflag=1;
}
}
(3)实验结果(截屏图):
(4).程序源代码(Matlab语言):
functiony=Newton(f,df,x0,eps,M)
d=0;
fork=1:
M
iffeval(df,x0)==0
d=2;break
else
x1=x0-feval(f,x0)/feval(df,x0);
end
e=abs(x1-x0);x0=x1;
ife<=eps&abs(feval(f,x1))<=eps
d=1;break
end
end
ifd==1
y=x1;
elseifd==0
y='迭代失败'
else
y='奇异'
end
functiony=df(x)
y=7*power(x,6)-112*power(x,3);
end
functiony=f(x)
y=power(x,7)-28*power(x,4)+14;
End
4.1计算结果打印:
(5)问题讨论与分析:
Newton迭代法是一个二阶收敛迭代式,他的几何意义
是
的切线与x轴的交点,故也称为切线法。
它是平方收敛的,但它是局部收敛的,即要求初始值与方程的根充分接近,所以在计算过程中需要先确定初始值。
本题在讨论中,讨论了区间(0.1,1.9)两端点是否能作为Newton迭代的初值,结果发现0.1不满足条件,而1.9满足,能作为初值
。
二,已知函数值如下表:
x
1
2
3
4
5
f(x)
0
0.693147
1.0986123
1.3862944
1.6094378
x
6
7
8
9
10
f(x)
1.7917595
1.9459101
2.079445
2.1972246
2.3025851
f'(x)
f'
(1)=1
f'(10)=0.1
试用三次样条插值求f(4.563)和f'(4.563)的近似值。
解:
(1)解题方法及理论依据:
公式1
公式2
三弯矩方程:
其矩阵形式:
其中
对于上述矩阵方程,由于其系数矩阵为三对角矩阵,所以用追赶法来解。
解三对角方程组的追赶法:
设线性方程组Ax=b的系数矩阵A为三对角矩阵
这时A可以作如下LU分解:
(1)LU分解:
首先
,对i=2,3,...,n,计算
(2)解Lz=b:
首先
对i=2,3,...,n,计算
(3)解Ux=z:
首先
对i=n-1,n-2,....,2,1,计算
最后把追赶法求得的唯一解
代入上面的公式1和公式2,即可求得
。
(2)计算机程序(C语言):
#include
#include
voidmain()
{
floatX[10]={1,2,3,4,5,6,7,8,9,10};//输入节点
floatf[10]={0,0.69314718,1.0986123,1.3862944,1.6094378,1.7917595,1.94594101,2.079445,2.1972246,2.3025851};//输入节点处的函数值
floath[10];
floatu[9],v[9],c[9];
floatr[10],z[10],l[10],a[10],d[10],m[10];
inti;
floatx=4.563;
floatS,s;
for(i=0;i<=8;i++)
h[i]=X[i+1]-X[i];//计算步长h
d[0]=6*((f[0]-f[1])/(X[0]-X[1])-1)/h[1];
d[8]=6*(0.1-(f[8]-f[9])/(X[8]-X[9]))/h[8];
u[8]=1;
v[0]=1;
for(i=0;i<=7;i++)
u[i]=h[i]/(h[i]+h[i+1]);//计算三弯矩方程组中的系数u
for(i=1;i<=8;i++)
{
v[i]=1-u[i];
d[i]=6*((f[i+1]-f[i])/(X[i+1]-X[i])-(f[i]-f[i-1])/(X[i]-X[i-1]))/(h[i-1]+h[i]);//计算三弯矩方程组中右边的常数项d
}
r[0]=2;//从此处开始用追赶法解三对角方程组
for(i=1;i<=9;i++)
{
l[i]=u[i]/r[i-1];
r[i]=2-l[i]*v[i-1];
}
z[0]=d[0];
for(i=1;i<=9;i++)
z[i]=d[i]-l[i]*z[i-1];
m[9]=z[9]/2;
for(i=8;i>=0;i--)
m[i]=(z[i]-v[i]*m[i+1])/2;//此处为前面用追赶法解方程组求得的结果
i=(int)x-1;
S=pow(X[i+1]-x,3)/(6*h[i+1])*m[i]+pow(x-X[i],3)/(6*h[i+1])*m[i+1]+(f[i]-(m[i]*pow(h[i+1],2))/6)*(X[i+1]-x)/h[i+1]+(f[i+1]-(m[i+1]*pow(h[i+1],2)/6))*(x-X[i])/h[i+1];//计算三次样条插值函数s(x)
s=-pow(X[i+1]-x,2)*m[i]/(2*h[i+1])+pow(x-X[i],2)*m[i+1]/(2*h[i+1])-(f[i]-m[i]*pow(h[i+1],2)/6)/h[i+1]+(f[i+1]-m[i+1]*pow(h[i+1],2)/6)/h[i+1];//计算s(x)的导数
printf("f(4.536)=%f\n",S);//输出f(4.536)的结果
printf("f'(4.536)=%f\n",s);//输出f'(4.536)的结果
}
(3)实验结果(截屏图):
(4).程序源代码(Matlab语言):
functionQ=san(ssss,p)
Q=zeros(2,1);
x=[1;2;3;4;5;6;7;8;9;10];
y=[0;0.69314718;1.0986123;1.3862944;1.6094378;1.7917595;1.9459101;2.079445;2.1972246;2.3025851];
h=zeros(10,1);
d=zeros(10,1);
u=zeros(10,1);
v=zeros(10,1);
r=zeros(10,1);
l=zeros(10,1);
z=zeros(10,1);
m=zeros(10,1);
fort=1:
1:
9;
h(t)=x(t+1)-x(t);
end
d
(1)=6/h
(1)*((y
(2)-y
(1))/h
(1)-1);
d(10)=6/h(9)*(0.1-(y(10)-y(9))/h(9));
fort=1:
1:
8
u(t+1)=h(t)/(h(t)+h(t+1));
v(t+1)=1-u(t+1);
d(t+1)=6/(h(t)+h(t+1))*((y(t+2)-y(t+1))/(x(t+2)-x(t+1))-(y(t+1)-y(t))/(x(t+1)-x(t)));
end
u(10)=1;v
(1)=1;r
(1)=d
(1);
fort=2:
1:
10
l(t)=u(t)/r(t-1);
r(t)=d(t)-l(t)*v(t-1);
end
z
(1)=d
(1);
fort=2:
1:
10
z(t)=d(t)-l(t)*z(t-1);
end
m(10)=z(10)/r(10);
fort=9:
-1:
1
m(t)=(z(t)-v(t)*m(t+1))/r(t);
end
fort=1:
1:
10
ifp>=t&&p<(t+1)
Q(:
1)=feval(ssss,p,t,x,m,h,y);break
end
end
functionQ=ssss(p,t,x,m,h,y)
Q=zeros(2,1);
Q(1,1)=((power((x(t+1)-p),3)*m(t)+power((p-x(t)),3)*m(t+1))/6+(y(t)-m(t)*h(t)*h(t)/6)*(x(t+1)-p)+(y(t+1)-m(t+1)*h(t)*h(t)/6)*(p-x(t)))/h(t);
Q(2,1)=(-(power((x(t+1)-p),2)*m(t)+power((p-x(t)),2)*m(t+1))/2+(y(t)-m(t)*h(t)*h(t)/6)+(y(t+1)-m(t+1)*h(t)*h(t)/6))/h(t);
End
4.1计算结果打印:
(5)问题讨论与分析:
样条插值效果比Lagrange插值好,近似误差较小,光滑性较好。
本题的对任意划分的三弯矩插值法可以解决非等距节点的一般性问题。
其实编程时对于数组的大小可以弄大点,没必要仔细考虑各数组变量的精确元素个数。
同时,编程时也不需要求程序中计算式与原基本公式下标完全一致(例如本题中用追赶法解三对角方程组时计算式中的变量名和书上就不一致),只须数组元素的对应位置(即带入地址)一致即可。
三,用Romberg算法求
(允许误差ε=0.00001)。
解:
(1)解题方法及理论依据:
所谓Romberg算法就是采用公式
、
、
、
、
按照如下图所示的加工流程计算高精度的romberg积分值
。
逐行计算,算完前五行后得
,若
(准许误差),就把
作为积分近似值,否则再算第六行得
,若
就把
作为积分近似值,否则.....,直到
时停止计算,并把
作为积分近似值。
(2)计算机程序(C语言):
#include
#include
voidmain()
{
floatf(floatx);//函数声明
floata=1,b=3,h,m;
intn,k;
floatR[10],S[10],T[10],C[10];
for(n=1;n<=10;n++)
{m=0;
h=(b-a)/(pow(2,(n-1)));//求步长h
for(k=0;k<(pow(2,(n-1)));k++)
m=m+f(a+(k+0.5)*h);
T[0]=(b-a)*(f(a)+f(b))/2;
T[n]=0.5*(T[n-1]+h*m);//求T
S[n-1]=(4*T[n]-T[n-1])/3;//求s
C[n-2]=(16*S[n-1]-S[n-2])/15;//求C
R[n-3]=(64*C[n-2]-C[n-3])/63;//求R
if(fabs(R[n-3]-R[n-4])<0.00001)//误差控制
break;
}
printf("%f\n",R[n-3]);//输出最终计算结果
}
floatf(floatx)//函数定义
{
floatf;
f=pow(3,x)*pow(x,1.4)*(5*x+7)*sin(x*x);
return(f);
}
(3)实验结果(截屏图):
(4)程序源代码(Matlab语言):
function[T,n]=mromb(f,a,b,eps)
ifnargin<4,eps=1e-6;end
h=b-a;
R(1,1)=(h/2)*(feval(f,a)+feval(f,b));
n=1;J=0;err=1;
while(err>eps)
J=J+1;h=h/2;S=0;
fori=1:
n
x=a+h*(2*i-1);
S=S+feval(f,x);
end
R(J+1,1)=R(J,1)/2+h*S;
fork=1:
J
R(J+1,k+1)=(4^k*R(J+1,k)-R(J,k))/(4^k-1);
end
err=abs(R(J+1,J+1)-R(J+1,J));
n=2*n;
end
R;
T=R(J+1,J+1);
4.1计算结果打印:
(5)问题讨论与分析:
Romberg算法作为一种求积方法加速的这么一种方法,在复化梯形公式的基础上,精确度提高了并且收敛速度加快了。
本题中求公式
中的
这一部分,其实也可以像被积函数那样做成一个函数来实现。
Romberg算法的缺点是:
对函数的光滑性要求较高,在计算新分点的值时,这些数值的个数成倍增加。
四,用定步长四阶Runge-Kutta法求解:
h=0.0005,打印yi(0.025),yi(0.045),yi(0.085),yi(0.1),(i=1,2,3)
解:
(1)解题方法及理论依据:
先不按Taylor公式展开,而是先写成tn处附近的值的线性组合(有待定常数)再按Taylor公式展开,然后确定待定常数,这就是Runge-Kutta法的思想方法。
定步长四阶Runge-Kutta公式的古典形式:
(2)计算机程序(C语言):
#include
#include
floatf1(floatx,floaty,floatz)//函数定义
{floatf;
f=1;
return(f);
}
floatf2(floatx,floaty,floatz)//函数定义
{floatf;
f=z;
return(f);
}
floatf3(floatx,floaty,floatz)//函数定义
{floatf;
f=1000-1000*y-100*z;
return(f);
}
voidmain()
{floath=0.0005,k1,k2,k3,k4,n1,n2,n3,n4,m1,m2,m3,m4,x1=0,x2,y1=0,y2,z1=0,z2,t=0;
inti=0;
//Runge-Kutta算法
do
{t=h*i;
if(i>=1)x1=x2,y1=y2,z1=z2;
//四阶古典Runge-Kutta公式
k1=h*f1(x1,y1,z1),
m1=h*f2(x1,y1,z1),
n1=h*f3(x1,y1,z1),
k2=h*f1(x1+k1/2,y1+m1/2,z1+n1/2),
m2=h*f2(x1+k1/2,y1+m1/2,z1+n1/2),
n2=h*f3(x1+k1/2,y1+m1/2,z1+n1/2),
k3=h*f1(x1+k2/2,y1+m2/2,z1+n2/2),
m3=h*f2(x1+k2/2,y1+m2/2,z1+n2/2),
n3=h*f3(x1+k2/2,y1+m2/2,z1+n2/2),
k4=h*f1(x1+k3,y1+m3,z1+n3),
m4=h*f2(x1+k3,y1+m3,z1+n3),
n4=h*f3(x1+k3,y1+m3,z1+n3),
x2=x1+(k1+2*k2+2*k3+k4)/6,
y2=y1+(m1+2*m2+2*m3+m4)/6,
z2=z1+(n1+2*n2+2*n3+n4)/6;
//输出结果
if(i==49)
printf("y1(%g)=%fy2(%g)=%fy3(%g)=%f\n",h*(i+1),x2,h*(i+1),y2,h*(i+1),z2);
if(i==89)
printf("y1(%g)=%fy2(%g)=%fy3(%g)=%f\n",h*(i+1),x2,h*(i+1),y2,h*(i+1),z2);
if(i==169)
printf("y1(%g)=%fy2(%g)=%fy3(%g)=%f\n",h*(i+1),x2,h*(i+1),y2,h*(i+1),z2);
if(i==199)
printf("y1(%g)=%fy2(%g)=%fy3(%g)=%f\n",h*(i+1),x2,h*(i+1),y2,h*(i+1),z2);
i++;
}
while(i<=200);
}
(3)实验结果(截屏图):
(4)程序源代码(Matlab语言):
functionY=R_K(df1,a,b,h)
m=(b-a)/h;
Y=zeros(3,1);
S=zeros(3,1);
K=zeros(3,4);
x=a;y1=a;y2=a;y3=a;
forn=1:
m
K(:
1)=feval(df1,x,y1,y2,y3);
x=x+0.5*h;S(:
1)=Y+0.5*h.*K(:
1);
y1=S(1,1);y2=S(2,1);y3=S(3,1);
K(:
2)=feval(df1,x,y1,y2,y3);
S(:
1)=Y+0.5*h.*K(:
2);
y1=S(1,1);y2=S(2,1);y3=S(3,1);
K(:
3)=feval(df1,x,y1,y2,y3);
x=x+0.5*h;S(:
1)=Y+h.*K(:
3);
y1=S(1,1);y2=S(2,1);y3=S(3,1);
K(:
4)=feval(df1,x,y1,y2,y3);
Y=Y+h.*(K(:
1)+2.*K(:
2)+2.*K(:
3)+K(:
4))/6;
end
functionZ=df1(x,y1,y2,y3)
Z=zeros(3,1);
Z
(1)=0*x+0*y1+0*y2+0*y3+1;
Z
(2)=0*x+0*y1+0*y2+1*y3;
Z(3)=0*x+0*y1-1000*y2-100*y3+1000;
End
4.1计算结果打印:
(5)问题讨论与分析:
Runge-Kutta方法的优点:
精度高,不必用别的方法求开始几点的函数值。
程序简单,存储量少。
方法稳定。
五,已知A与b
用列主元消去法求解Ax=b。
解:
(1)解题方法及理论依据:
列主元素消去法是建立在Gauss消去法基础上的,Gauss消去法是在较强条件下进行的,它要求主元
另外,当主元的绝对值很小时,会损失精度。
为了避免Gauss消去法的数值不稳定性,在每一步消元之前增加一个选主元的过程,我们选取绝对值尽量大的元素作为主元,在每一次消元之前,就按列选择一次主元。
完成了n-1次消元之后,再回代的方法称为Gauss列主元素消元法。
(2)计算机程序(C语言):
#include
#include
#defineN9
doubleLZXY(doublea[N][N+1],doublex[N])//用列主元素消去法求解线性方程组的解的函数
{
inti,j,k,r;
doublec,l,t[N];
for(k=0;k { r=k; for(i=k;i if(fabs(a[i][k])>fabs(a[r][k]))r=i;//选绝对值最大的元素作为主元 if(r! =k) { for(j=k;j<=N;j++)//对调r,k两行 { c=a[k][j]; a[k][j]=a[r][j]; a[r][j]=c; } } for(i=k+1;i { l=a[i][k]/a[k][k]; for(j=k;j<=N;j++) a[i][j]=a[i][j]-l*a[k][j]; } } x[N-1]=a[N-1][N]/a[N-1][N-1];//回代求解 for(i=N-2;i>=0;i--) { t[i]=a[i][N]; for(j=i+1;j<=N-1;j++) t[i]=t[i]-a[i][j]*x[j]; x[i]=t[i]/a[i][i]; } return (1); } voidmain() { inti; doubledet,x[N]; doublea[N][N+1]={ {12.38412,2.115237,-1.061074,1.112336,-0.113584,0.718719,1.742382,3.067813,-2.031743,2.1874369}, {2.115237,19.141823,-3.125432,-1.012345,2.189736,1.563849,-0.784165,1.112348,3.123124,33.992318}, {-1.061074,-3.125432,15.567914,3.123848,2.031454,1.836742,-1.056781,0.336993,-1.010103,-25.
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 13 级数 分析 实验 报告