cout<<"是否输入其余要求x的值[是
(1),否(0)]:
";
cin>>q;
}
system("cls");
}
voidmain()
{
lagrange();
}
实验结果:
实验二
实验题目:
用不同的方法计算积分
取不同的步长h,分别用复合梯形公式及复合辛普森公式计算积分,给出误差中关于h的函数,并与积分精确值比较两个公式的精度,是否存在一个最小的h,使得精度不能再被改善。
程序代码:
复合梯形公式:
functionI=T_quad(x,y)
n=length(x);m=length(y);
ifn~=m
error('ThelengthsofXandYmustbeequal');
return;
end
h=(x(n)-x
(1))/(n-1);a=[12*ones(1,n-2)1];
I=h/2*sum(a.*y);
复合辛普森公式:
functionI=S_quad(x,y)
n=length(x);m=length(y);
ifn~=m
error('ThelengthsofXandYmustbeequal');
return;
end
ifrem(n-1,2)~=0
I=T_quad(x,y);
return;
end
N=(n-1)/2;h=(x(n)-x
(1))/N;a=zeros(1,n);
fork=1:
N
a(2*k-1)=a(2*k-1)+1;a(2*k)=a(2*k)+4;
a(2*k+1)=a(2*k+1)+1;
end
I=h/6*sum(a.*y);
测试数据:
复合梯形测试数据:
>>x=0.00001:
0.0001:
0.99999
>>y=sqrt(x).*log(x)
>>I=T_quad(x,y)
复合辛普森测试数据:
>>x=0.00001:
0.0001:
0.99999
>>y=sqrt(x).*log(x)
>>I=S_quad(x,y)
实验结果:
复合梯形实验结果:
>>x=0.00001:
0.0001:
0.99999
x=
Columns1through8
0.00000.00010.00020.00030.00040.00050.00060.0007
Columns9through16
…………
Columns9993through10000
0.99920.99930.99940.99950.99960.99970.99980.9999
>>y=sqrt(x).*log(x)
y=
Columns1through8
-0.0364-0.0956-0.1227-0.1422-0.1579-0.1712-0.1828-0.1932
…………
Columns9993through10000
-0.0008-0.0007-0.0006-0.0005-0.0004-0.0003-0.0002-0.0001
>>I=T_quad(x,y)
I=
-0.4444
I=
-0.4444
复合辛普森实验结果:
>>x=0.00001:
0.0001:
0.99999
x=
Columns1through8
0.00000.00010.00020.00030.00040.00050.00060.0007
…………
Columns9993through10000
0.99920.99930.99940.99950.99960.99970.99980.9999
>>y=sqrt(x).*log(x)
y=
Columns1through8
-0.0364-0.0956-0.1227-0.1422-0.1579-0.1712-0.1828-0.1932
…………
Columns9993through10000
-0.0008-0.0007-0.0006-0.0005-0.0004-0.0003-0.0002-0.0001
>>I=S_quad(x,y)
I=
-0.4444
I=
-0.4444
实验三
实验题目:
用LU分解和列主元消去法解线性方程组
输出Ax=b中系数A=LU分解的矩阵L和U,解向量x及detA;列主元法的行交换次序,解向量x及detA;比较两种方法所得的结果。
程序代码:
//解线性方程组
#include
#include
#include
//----------------------------------------------全局变量定义区
constintNumber=15;//方程最大个数
doublea[Number][Number],b[Number],copy_a[Number][Number],copy_b[Number];//系数行列式
intA_y[Number];//a[][]中随着横坐标增加列坐标的排列顺序,如a[0][0],a[1][2],a[2][1]...则A_y[]={0,2,1...};
intlenth,copy_lenth;//方程的个数
doublea_sum;//计算行列式的值
char*x;//未知量a,b,c的载体
//----------------------------------------------函数声明区
voidinput();//输入方程组
voidprint_menu();//打印主菜单
intchoose();//输入选择
voidcramer();//Cramer算法解方程组
voidgauss_row();//Gauss列主元解方程组
voidDoolittle();//用Doolittle算法解方程组
intDoolittle_check(doublea[][Number],doubleb[Number]);//判断是否行列式>0,若是,调整为顺序主子式全>0
voidxiaoqu_u_l();//将行列式Doolittle分解
voidcalculate_u_l();//计算Doolittle结果
double&calculate_A(intn,intm);//计算行列式
doublequanpailie_A();//根据列坐标的排列计算的值,如A_y[]={0,2,1},得sum=a[0][A_y[0]]*a[1][A_y[1]]*a[2][A_y[2]]=a[0][0]*a[1][2]*a[2][1];
voidexchange(intm,inti);//交换A_y[m],A_y[i]
voidexchange_lie(intj);//交换a[][j]与b[];
voidexchange_hang(intm,intn);//分别交换a[][]和b[]中的m与n两行
voidgauss_row_xiaoqu();//Gauss列主元消去法
voidgauss_calculate();//根据Gauss消去法结果计算未知量的值
voidexchange_a_lie(intm,intn);//交换a[][]中的m和n列
voidexchange_x(intm,intn);//交换x[]中的x[m]和x[n]
voidrecovery();//恢复数据
//主函数
voidmain()
{
intflag=1;
input();//输入方程
while(flag)
{
print_menu();//打印主菜单
flag=choose();//选择解答方式
}
}
//函数定义区
voidprint_menu()
{
system("cls");
cout<<"------------方程系数和常数矩阵表示如下:
\n";
for(intj=0;jcout<<"系数"<cout<<"\t常数";
cout<for(inti=0;i{
for(j=0;jcout<:
left)<cout<<"\t"<
}
cout<<"-----------请选择方程解答的方案----------";
cout<<"\n1.Gauss列主元消去法";
cout<<"\n2.Doolittle分解法";
cout<<"\n3.退出";
cout<<"\n输入你的选择:
";
}
voidinput()
{inti,j;
cout<<"方程的个数:
";
cin>>lenth;
if(lenth>Number)
{
cout<<"Itistoobig.\n";
return;
}
x=newchar[lenth];
for(i=0;ix[i]='a'+i;
//输入方程矩阵
//提示如何输入
cout<<"====================================================\n";
cout<<"请在每个方程里输入"<\n";
cout<<"例:
\n方程:
a";
for(i=1;i{
cout<<"+"<
}
cout<<"=10\n";
cout<<"应输入:
";
for(i=0;icout<
cout<<"10\n";
cout<<"==============================\n";
//输入每个方程
for(i=0;i{
cout<<"输入方程"<
";
for(j=0;jcin>>a[i][j];
cin>>b[i];
}
//备份数据
for(i=0;ifor(j=0;jcopy_a[i][j]=a[i][j];
for(i=0;icopy_b[i]=b[i];
copy_lenth=lenth;
}
//输入选择
intchoose()
{
intchoice;charch;
cin>>choice;
switch(choice)
{
case1:
cramer();break;
case2:
gauss_row();break;
case3:
guass_all();break;
case4:
Doolittle();break;
case5:
return0;
default:
cout<<"输入错误,请重新输入:
";
choose();
break;
}
cout<<"\n是否换种方法求解(Y/N):
";
cin>>ch;
if(ch=='n'||ch=='N')return0;
recovery();
cout<<"\n\n\n";
return1;
}
//高斯列主元排列求解方程
voidgauss_row()
{
inti,j;
gauss_row_xiaoqu();//用高斯列主元消区法将系数矩阵变成一个上三角矩阵
for(i=0;i{
for(j=0;jcout<cout<}
if(a[lenth-1][lenth-1]!
=0)
{
cout<<"系数行列式不为零,方程有唯一的解:
\n";
gauss_calculate();
for(i=0;i{
cout<}
}
else
cout<<"系数行列式等于零,方程没有唯一的解.\n";
}
voidgauss_row_xiaoqu()//高斯列主元消去法
{
inti,j,k,maxi;doublelik;
cout<<"用Gauss列主元消去法结果如下:
\n";
for(k=0;k{
j=k;
for(maxi=i=k;iif(a[i][j]>a[maxi][j])maxi=i;
if(maxi!
=k)
exchange_hang(k,maxi);//
for(i=k+1;i{
lik=a[i][k]/a[k][k];
for(j=k;ja[i][j]=a[i][j]-a[k][j]*lik;
b[i]=b[i]-b[k]*lik;
}
}
}
voidDoolittle()//Doolittle消去法计算方程组
{
doubletemp_a[Number][Number],temp_b[Number];inti,j,flag;
for(i=0;ifor(j=0;jtemp_a[i][j]=a[i][j];
flag=Doolittle_check(temp_a,temp_b);
if(flag==0)cout<<"\n行列式为零.无法用Doolittle求解.";
xiaoqu_u_l();
calculate_u_l();
cout<<"用Doolittle方法求得结果如下:
\n";
for(i=0;i{
for(j=0;x[j]!
='a'+i&&jcout<}
}
voidcalculate_u_l()//计算Doolittle结果
{inti,j;doublesum_ax=0;
for(i=0;i{
for(j=0,sum_ax=0;j
sum_ax+=a[i][j]*b[j];
b[i]=b[i]-sum_ax;
}
for(i=lenth-1;i>=0;i--)
{
for(j=i+1,sum_ax=0;jsum_ax+=a[i][j]*b[j];
b[i]=(b[i]-sum_ax)/a[i][i];
}
}
voidxiaoqu_u_l()//将行列式按Doolittle分解
{inti,j,n,k;doubletemp;
for(i=1,j=0;ia[i][j]=a[i][j]/a[0][0];
for(n=1;n{//求第n+1层的上三角矩阵部分即U
for(j=n;j{for(k=0,temp=0;ktemp+=a[n][k]*a[k][j];
a[n][j]-=temp;
}
for(i=n+1;i{for(k=0,temp=0;ktemp+=a[i][k]*a[k][n];
a[i][n]=(a[i][n]-temp)/a[n][n];
}
}
}
intDoolittle_check(doubletemp_a[][Number],doubletemp_b[Number])//若行列式不为零,将系数矩阵调整为顺序主子式大于零
{
inti,j,k,maxi;doublelik,temp;
for(k=0;k{
j=k;
for(maxi=i=k;iif(temp_a[i][j]>temp_a[maxi][j])maxi=i;
if(maxi!
=k)
{exchange_hang(k,maxi);
for(j=0;j{temp=temp_a[k][j];
temp_a[k][j]=temp_a[maxi][j];
temp_a[maxi][j]=temp;
}
}
for(i=k+1;i{
lik=temp_a[i][k]/temp_a[k][k];
for(j=k;jtemp_a[i][j]=temp_a[i][j]-temp_a[k][j]*lik;
temp_b[i]=temp_b[i]-temp_b[k]*lik;
}
}
if(temp_a[lenth-1][lenth-1]==0)return0;
return1;
}
voidexchange_hang(intm,intn)//交换a[][]中和b[]两行
{
intj;doubletemp;
for(j=0;j{temp=a[m][j];
a[m][j]=a[n][j];
a[n][j]=temp;
}
temp=b[m];
b[m]=b[n];
b[n]=temp;
}
voidexchange(intm,inti)//交换A_y[m],A_y[i]
{inttemp;
temp=A_y[m];
A_y[m]=A_y[i];
A_y[i]=temp;
}
voidexchange_lie(intj)//交换未知量b[]和第i列
{doubletemp;inti;
for(i=0;i{temp=a[i][j];
a[i][j]=b[i];
b[i]=temp;
}
}
voidexchange_a_lie(intm,intn)//交换a[]中的两列
{doubletemp;inti;
for(i=0;i{temp=a[i][m];
a[i][m]=a[i][n];
a[i][n]=temp;
}
}
voidexchange_x(intm,intn)//交换未知量x[m]与x[n]
{chartemp;
temp=x[m];
x[m]=x[n];
x[n]=temp;
}
voidrecovery()//用其中一种方法求解后恢复数据以便用其他方法求解
{
for(inti=0;ifor(intj=0;ja[i]