计算方法作业.docx
《计算方法作业.docx》由会员分享,可在线阅读,更多相关《计算方法作业.docx(14页珍藏版)》请在冰豆网上搜索。
![计算方法作业.docx](https://file1.bdocx.com/fileroot1/2023-6/8/88784520-5284-4a94-810c-06e6b777f4c2/88784520-5284-4a94-810c-06e6b777f4c21.gif)
计算方法作业
计算方法练习
第一题、
梯形公式:
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);
}
运行结果: