计算方法作业.docx
《计算方法作业.docx》由会员分享，可在线阅读，更多相关《计算方法作业.docx（14页珍藏版）》请在冰豆网上搜索。
计算方法作业
计算方法练习
第一题、
梯形公式：
C程序：
#include
#include
#definee2.718281828
floatfun（floatx）{
floatf;
f=（1-pow（e,（-x）））/x;
returnf;
}
voidmain（）{
inti,n;
floata,b;
floath,m,T,s=0;
printf（"pleaseinputa="）;
scanf（"%f",&a）;
printf（"pleaseinputb="）;
scanf（"%f",&b）;
printf（"pleaseinputn="）;
scanf（"%d",&n）;
if（b>a）{
h=（b-a）/n;
m=（b-a）/n;
}else{
printf（"errorinput!
"）;
}
for（i=1;i<=n-1;i++）{
m=m*i;
s=s+2*fun（a+m）;
}
T=h/2*（fun（a）+s+fun（b））;
printf（"T=%f\n",T）;
}
运行结果：
辛普森公式：
#include
#include
#definee2.718281828
floatfun（floatx）{
floatf;
f=（1-pow（e,（-x）））/x;
returnf;
}
voidmain（）{
inti,n;
floata,b,h,m,sum=0,S;
printf（"pleaseinputa="）;
scanf（"%f",&a）;
printf（"pleaseinputb="）;
scanf（"%f",&b）;
printf（"pleaseinputn="）;
scanf（"%d",&n）;
if（b>a）{
h=（b-a）/n;
m=h/2;
}else{
printf（"errorinput!
"）;
}
for（i=1;i<=n-1;i++）{
sum=sum+4*fun（a+m*（2*i-1））+2*fun（a+m*（2*i））;
}
S=h/6*（fun（a）+sum+fun（b））;
printf（"S=%f\n",S）;
}
运行结果：
（2）
梯形公式：
C程序：
#include
#include
floatfun（floatx）{
floatf;
f=x/（4+x*x）;
returnf;
}
voidmain（）{
inti,n;
floata,b;
floath,m,T,s=0;
printf（"pleaseinputa="）;
scanf（"%f",&a）;
printf（"pleaseinputb="）;
scanf（"%f",&b）;
printf（"pleaseinputn="）;
scanf（"%d",&n）;
if（b>a）{
h=（b-a）/n;
m=（b-a）/n;
}else{
printf（"errorinput!
"）;
}
for（i=1;i<=n-1;i++）{
m=m*i;
s=s+2*fun（a+m）;
}
printf（"b-a=%f\n",b-a）;
printf（"h=%f\n",h）;
printf（"m=%f\n",m）;
printf（"s=%f\n",s）;
T=h/2*（fun（a）+s+fun（b））;
printf（"T=%f\n",T）;
}
运行结果：
辛普森公式：
#include
#include
floatfun（floatx）{
floatf;
f=x/（4+x*x）;
returnf;
}
voidmain（）{
inti,n;
floata,b,h,m,sum=0,S;
printf（"pleaseinputa="）;
scanf（"%f",&a）;
printf（"pleaseinputb="）;
scanf（"%f",&b）;
printf（"pleaseinputn="）;
scanf（"%d",&n）;
if（b>a）{
h=（b-a）/n;
m=h/2;
}else{
printf（"errorinput!
"）;
}
for（i=1;i<=n-1;i++）{
sum=sum+4*fun（a+m*（2*i-1））+2*fun（a+m*（2*i））;
}
S=h/6*（fun（a）+sum+fun（b））;
printf（"S=%f\n",S）;
}
运行结果：
第二题、
建立M文件guass1.m
functions=guassl（a,b,n）
h=（b-a）/n;
s=0.0;
form=0:
（1*n/2-1）
s=s+h*（guassf（a+h*（（1-1/sqrt（3））+2*m））+guassf（a+h*（（1+1/sqrt（3））+2*m）））;
end
s;
I=int（'sin（x）',0,1）;
c=（I-s）/I;d=vpa（c,10）;
建立M文件guassf.m
functiony=guassf（x）
y=sqrt（x）*exp（x）;
建立M文件guassf.m
functiony=guassf（x）
y=sqrt（-x*x）*exp（1+x*x）;
运行：
s=guassl（0,inf,4）
s=guassl（0,inf,6）
s=guassl（-inf,inf,4）
s=guassl（-inf,inf,6）
第三题、
建立M文件my_LU.m
functionmy_LU（a,b）
n=length（a）;
l=zeros（n,n）;
u=zeros（n,n）;
fori=1:
n
l（i,i）=1;
end
u（1,1:
n）=a（1,1:
n）;
l（2:
n,1）=a（2:
n,1）./u（1,1）;
forr=2:
n
fori=r:
n
u（r,i）=a（r,i）-sum（l（r,1:
r-1）.*（u（1:
r-1,i））'）;
end
fori=r+1:
n
if（r~=n）
l（i,r）=（a（i,r）-sum（l（i,1:
r-1）.*（u（1:
r-1,r）））'）./u（r,r）;
end
end
end
L=l,U=u
y
（1）=b
（1）;
fori=2:
n
y（i）=b（i）-sum（l（i,1:
i-1）.*y（1:
i-1））;
end
x（n）=y（n）/u（n,n）;
fori=n-1:
-1:
1
x（i）=（y（i）-sum（u（i,i+1:
n）.*x（i+1:
n）））./u（i,i）;
end
x=x'
end
在窗口运行：
>>a=[234;119;12-6]
a=
234
119
12-6
>>b=[021]
b=
021
>>my_LU（a,b）
L=
1.000000
0.50001.00000
0.5000-1.00001.0000
U=
2.00003.00004.0000
0-0.50007.0000
00-1.0000
x=
75
-46
-3
第四题、
MATLAB代码如下：
>>A=[3,2,1;
2,1,0;
1,0,3;];
b=[5,3,4];
n=length（b）;%方程个数n
G=zeros（n,n）;
G（1,1）=sqrt（A（1,1））;
G（2:
n,1）=A（2:
n,1）/G（1,1）;
forj=2:
n-1
G（j,j）=sqrt（A（j,j）-sum（G（j,1:
j-1）.^2））;
fori=j+1:
n
G（i,j）=（A（i,j）-sum（G（i,1:
j-1）.*G（j,1:
j-1）））/G（j,j）;
end
end
G（n,n）=sqrt（A（n,n）-sum（（G（n,1:
n-1））.^2））;
G
%--------用A=GG’分解求解方程组Ax=b（平方根法）----
%Gy=b，G'x=y;
x=zeros（n,1）;%未知向量
y=zeros（n,1）;%中间向量
y
（1）=b
（1）/G（1,1）;
fori=2:
n
y（i）=（b（i）-sum（G（i,1:
i-1）'.*y（1:
i-1）））/G（i,i）;
end
y
%---------由G'x=y求出x----------
%方法类似于U由x=y解出x
G=G';
x（n）=y（n）/G（n,n）;
fori=n-1:
-1:
1
x（i）=（y（i）-sum（G（i,i+1:
n）'.*x（i+1）））/G（i,i）;
end
x
G=
1.7321+0.0000i0.0000+0.0000i0.0000+0.0000i
1.1547+0.0000i0.0000+0.5774i0.0000+0.0000i
0.5774+0.0000i0.0000+1.1547i2.0000+0.0000i
y=
2.8868+0.0000i
0.0000+0.5774i
0.8333+0.0000i
x=
1.8333
-0.1667
0.4167
第五题、
（1）、牛顿法：
functionrtn=newton1（fx,dfx,x0,tol,N）
x=x0;f0=eval（fx）;df0=eval（dfx）;
n=0;
disp（'[nxnxn+1fn+1]'）;
whilen<=N
x1=x0-f0/df0;
x=x1;f1=eval（fx）;
X=[n,x0,x1,f1];
disp（X）;
ifabs（x0-x1）
fprintf（'Theprocedurewassuccessful.'）
return
else
n=n+1;
x0=x1;f0=f1;
end
end
ifn==N+1
fprintf（'themethodfailedafterNiterations.'）,
end
>>fx='x*x*x-3*x-1';
>>dfx='3*x*x-3';
newton1（fx,dfx,.1,10^（-5）,2）
[nxnxn+1fn+1]
00.1000-0.3374-0.0263
（弦截法）
C程序：
#include
#include
voidmain（）
{
floatf（floatx）;
floatroot（floatx1,floatx2）;
floatxpoint（floatx1,floatx2）;
floatx1,x2,f1,f2,x;
do
{
printf（"inputx1,x2:
\n"）;
scanf（"%f,%f",&x1,&x2）;
f1=f（x1）;
f2=f（x2）;
}while（f1*f2>=0）;
x=root（x1,x2）;
printf（"Arootofequationis%8.4f\n",x）;
}
floatroot（floatx1,floatx2）
{
floatxpoint（floatx1,floatx2）;
floatf（floatx）;
inti;
floatx,y,y1;
y1=f（x1）;
do
{
x=xpoint（x1,x2）;
y=f（x）;
if（y*y1>0）
{
y1=y;
x1=x;
}else
x2=x;
}while（fabs（y）>=0.0001）;
return（x）;
}
floatf（floatx）
{
floaty;
y=x*x*x-3*x-1;
return（y）;
}
floatxpoint（floatx1,floatx2）
{
floaty;
y=（x1*f（x2）-x2*f（x1））/（f（x2）-f（x1））;
return（y）;
}
运行结果：
升级会员 