数值分析编程题c语言汇总.docx
- 文档编号:8053216
- 上传时间:2023-01-28
- 格式:DOCX
- 页数:16
- 大小:24.08KB
数值分析编程题c语言汇总.docx
《数值分析编程题c语言汇总.docx》由会员分享,可在线阅读,更多相关《数值分析编程题c语言汇总.docx(16页珍藏版)》请在冰豆网上搜索。
数值分析编程题c语言汇总
数值分析实习报告
-1-
数值分析实习报告
上机实习题一一、题目:
b
与已知A12.38412,2.115237,-1.061074,1.112336,-0.113584,0.718719,1.742382,3.067813,-2.031743
2.115237,19.141823,-3.125432,-1.012345,2.189736,1.563849,-0.784165,1.112348,3.123124
-1.061074,-3.125432,15.567914,3.123848,2.031454,1.836742,-1.056781,0.336993,-1.010103
1.112336,-1.012345,3.123848,27.108437,4.101011,-3.741856,2.101023,-0.71828,-0.037585
A=-0.113584,2.189736,2.031454,4.101011,19.897918,0.431637,-3.111223,2.121314,1.784137
0.718719,1.563849,1.836742,-3.741856,0.431637,9.789365,-0.103458,-1.103456,0.238417
1.742382,-0.784165,-1.056781,2.101023,-3.111223,-0.103458,14.7138465,3.123789,-2.213474
3.067813,1.112348,0.336993,-0.71828,2.121314,-1.103456,3.123789,30.719334,4.446782
-2.031743,3.123124,-1.010103,-0.037585,1.784317,0.238417,-2.213474,4.446782,40.00001
b={2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4.719345,-5.6784392}
1.用Household变换,把A化为三对角阵B(并打印B)。
(0)=0,迭代9次)ω=1.4,X。
2.用超松弛法求解BX=b(取松弛因子3.用列主元素消去法求解BX=b。
二、解题方法的理论依据:
1、用Householder变换的理论依据
﹝1﹞令A0=A,a(ij)1=a(ij),已知Ar_1即Ar_1=a(ij)r
﹝2﹞Sr=sqrt(pow(a,2))
﹝3﹞a(r)=Sr*Sr+abs(a(r+1,r))*Sr
﹝4﹞y(r)=A(r_1)*u?
/a?
﹝5﹞Kr=(/2)*Ur的转置*Yr/a?
﹝6﹞Qr=Yr-Kr*Ur
﹝7﹞Ar=A(r-1)-(Qr*Ur的转置+Ur*Qr的转置)r=1,2,,……,n-2
2、用超松弛法求解
(m)(m-1)其基本思想:
在高斯方法已求出x,x的基础上,组合新的序列,从而加快收敛速度。
其算式:
B[i][i?
1]?
X[i?
1]?
B[i][i]?
X[i]?
B[i][i?
1]?
X[i?
1]?
b[i]
B[i][i?
1]B[i][i?
1]b[i]~?
X[i]?
?
X[i?
1]?
X0[i?
1]?
?
B[i][i]B[i][i]B[i][i]?
~X[i]?
w?
X[i]?
X0[i]?
?
X0[i]?
X[i]?
?
其中ω是超松弛因子,当ω>1时,可以加快收敛速度
3、用消去法求解
用追赶消去法求Bx=b的方法:
d1[i?
1]?
b[i]a1[i?
2]?
B[i?
1][i],
-2-
数值分析实习报告
]1i][i?
?
1[i1]?
B[ib1[?
1]?
B[i][i]c,,
q1[0]=0,u1[0]=0,
8?
,?
?
i]),i?
1,2ib1[]?
a1[i]?
q1[1q1[i]?
?
c[i](
u1[i}?
(d1[i]?
a1[i]?
u1[i])(b1[i]?
a1[i]?
q1[i?
1]),i?
1,2,?
?
?
9
x[9]=u1[9]
x[i]?
q1[i]?
x[i?
1]?
u1[i],i?
8,7,?
?
?
1
三、1.计算程序:
#includemath.h
#includestdio.h
#definege8
voidmain()
{
intsign(doublex);
doublea[][9]=
{
{12.38412,2.115237,-1.061074,1.112336,-0.113584,0.718719,1.742382,3.067813,-2.031743},
{2.115237,19.141823,-3.125432,-1.012345,2.189736,1.563849,-0.784165,1.112348,3.123124},
{-1.061074,-3.125432,15.567914,3.123848,2.031454,1.836742,-1.056781,0.336993,-1.010103},
{1.112336,-1.012345,3.123848,27.108437,4.101011,-3.741856,2.101023,-0.71828,-0.037585},
{-0.113584,2.189736,2.031454,4.101011,19.897918,0.431637,-3.111223,2.121314,1.784317},
{0.718719,1.563849,1.836742,-3.741856,0.431637,9.789365,-0.103458,-1.103456,0.238417},
{1.742382,-0.784165,-1.056781,2.101023,-3.111223,-0.103458,14.713846,3.123789,-2.213474},
{3.067813,1.112348,0.336993,-0.71828,2.121314,-1.103456,3.123789,30.719334,4.446782},
{-2.031743,3.123124,-1.010103,-0.037585,1.784317,0.238417,-2.213474,4.446782,40.00001}
};
doublek,h,s,w;
inti,j,n,m,g;
doubleu[9],x1[9],y[9],q[9],b1[9][10],x[9];
doubleb[9]=
{2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4.719345,-5.6784392};
for(j=0;j<7;++j)/*Household变换*/
{
s=0.0;
for(i=j+1;i<9;++i)
s=s+a[i][j]*a[i][j];
s=sqrt(s);
h=(a[j+1][j]>0)?
(s*s+s*a[j+1][j]):
(s*s-s*a[j+1][j]);
for(g=0;g<9;++g)
{
-3-
数值分析实习报告
if(g<=j)
u[g]=0;
elseif(g==j+1)
u[g]=a[j+1][j]+s*sign(a[j+1][j]);
elseu[g]=a[g][j];
}
for(m=0;m<9;++m)
{
y[m]=0;
for(n=0;n<9;++n)
y[m]=y[m]+a[m][n]*u[n];
y[m]=y[m]/h;
}
k=0;
for(i=0;i<9;++i)
k=k+u[i]*y[i];
k=0.5*k/h;
for(i=0;i<9;++i)
q[i]=y[i]-k*u[i];
for(n=0;n<9;++n)
for(m=0;m<9;++m)
a[m][n]=a[m][n]-(q[m]*u[n]+u[m]*q[n]);
}
printf(Household:
\n);
for(i=0;i<9;++i)
for(j=0;j<9;++j)
{
if(j%9==0)
printf(\
);
printf(%-9.5f,a[i][j]);
}
printf(\
);
w=1.4;/*超松弛法*/
for(i=0;i<9;i++)
x1[i]=0;
for(i=0;i<9;i++)
for(j=0;j<9;j++)
{
if(i==j)
b1[i][j]=0;
elseb1[i][j]=-a[i][j]/a[i][i];
}
for(i=0;i<9;i++)
b1[i][9]=b[i]/a[i][i];
-4-
数值分析实习报告
for(n=0;n<9;n++)
for(i=0;i<9;i++)
{
s=0;
for(j=0;j<9;j++)
s=s+b1[i][j]*x1[j];
s=s+b1[i][9];
x1[i]=x1[i]*(1-w)+w*s;
}
for(i=0;i<9;i++)
{
if(i==5)
printf(\
);
牰湩晴尨?
?
?
?
?
屦椬砬嬱嵩?
}
printf(\
);
u[0]=a[0][0];/*以下是消去法*/
y[0]=b[0];
for(i=1;i<9;i++)
{
q[i]=a[i][i-1]/u[i-1];
u[i]=a[i][i]-q[i]*a[i-1][i];
y[i]=b[i]-q[i]*y[i-1];
}
x[ge]=y[ge]/u[ge];
for(i=ge-1;i>=0;i--)
x[i]=(y[i]-a[i][i+1]*x[i+1])/u[i];
for(i=0;i<9;i++)
{
if(i==5)
printf(\
);
printf(x%d=%-10.6f,i,x[i]);
}
}
intsign(doublex)
{
intz;
z=(x>=(1e-6)?
1:
-1);
return(z);
}
2.打印结果:
Household:
12.38412-4.893080.000000.000000.000000.000000.000000.000000.00000
-5-
数值分析实习报告
-4.8930825.398426.494100.000000.000000.000000.000000.000000.00000
0.000006.4941020.611508.243930.000000.000000.000000.000000.00000
0.000000.000008.2439323.42284-13.880070.000000.000000.00000-0.00000
0.000000.000000.00000-13.8800729.698284.534500.000000.000000.00000
0.000000.000000.000000.000004.5345016.006124.881440.000000.00000
0.000000.000000.000000.000000.000004.8814426.01332-4.50363-0.00000
0.000000.000000.000000.000000.000000.00000-4.5036321.254064.50450
0.000000.000000.00000-0.000000.000000.00000-0.000004.5045014.53412
x0=1.073409x1=2.272579x2=-2.856601
x3=2.292514x4=2.112165x5=-6.422586
x6=1.357802x7=0.634259x8=-0.587042
x0=1.075799x1=2.275744x2=-2.855515
x3=2.293099x4=2.112634x5=-6.423838
x6=1.357923x7=0.634244x8=-0.587266
四、问题讨论:
此程序具有很好的通用性。
在GS方法的基础上,已经求出x的第m解,第m-1基础上,经过重新组合得新的序列,而在此新序列收敛速度加快。
上机实习题二
一、题目:
已知函数值如下表:
x
1
2
3
4
5
f(x)
0
0.6931478
1.0986123
1.3862944
1.6094378
x
6
7
8
9
10
f(x)
1.791795
1.9459101
2.079445
2.1972246
2.3025851
f'(x)
f'
(1)=1
f'(0)=0.1
试用三次样条插值求f(4.563)及f′(4.563)的近似值。
二、解题方法的理论依据:
任意划分的三弯矩插值法以及方程组解法中的三对角阵追赶算法。
应用三次样条插值法能够对函数产生很好的逼近效果。
而追赶算法又具有计算量少、方法简单、算法稳定的特点。
方法应用条件:
适用于求复杂函数在给定区间内某一点的函数值,给出函数f(x)在区间[a,b]中的n个插值点,并且给出函数在区间端点处的值。
三、1.计算程序:
#includestdio.h
#includemath.h
#definen11
-6-
数值分析实习报告
#definege11
voidmain()
{
inti,m;
doublee,s,E,q[12],u[12],y[12],c[12],w[12];
doubleb[12]={2,0,4.15888308,6.5916738,8.3177664,9.6566268,10.750557,11.6754606,
12.47667,13.1833476,13.8155106,14.0155106};
};
a[n][n-1]=4;
a[n][n]=2;
for(i=1;i<11;i++)
{
a[i][i-1]=1;
a[i][i]=4;
a[i][i+1]=1;
}
u[0]=a[0][0];/*消去法求c[i]*/
y[0]=b[0];
for(i=1;i<12;i++)
{
q[i]=a[i][i-1]/u[i-1];
u[i]=a[i][i]-q[i]*a[i-1][i];
y[i]=b[i]-q[i]*y[i-1];
}
c[ge]=y[ge]/u[ge];
for(i=ge-1;i>=0;i--)
c[i]=(y[i]-a[i][i+1]*c[i+1])/u[i];
牰湩晴尨请输入要插的值:
);
scanf(%lf,&E);
for(i=0;i<12;i++)
{
e=fabs(E-i);
if(e>=2)
w[i]=0;
elseif(e<=1)
w[i]=0.5*fabs(e*e*e)-e*e+2.0/3.0;
else
w[i]=(-1.0/6.0)*fabs(e*e*e)+e*e-2*fabs(e)+4.0/3.0;
}
s=0.0;
for(i=0;i<12;i++)
s=s+c[i]*w[i];
printf(
(%lf)=%lf,E,s);
printf(\
);
-7-
数值分析实习报告
牰湩晴尨请输入要求的导数的值:
);
scanf(%d,&m);
printf(
'(%d)=%lf\n,m,(c[m+1]-c[m-1])/2.0);
}
输出结果:
请输入要插的值:
4.563
f(4.563)=1.517932
请输入要求的导数的值:
4.563
f'(4.563)=0.249350
四、问题讨论:
在给均匀分划的插值函数x赋值时,由于使用for循环,误将x[i]=i+1写成x[i]=i,导致运算错误。
此程序具有一定通用性,对于任意划分的三弯矩插值法,只许改动x[i]即可。
求解方程组Mj时,要用到三对角方程组的追赶法(也称Thomas算法)。
变量较多,应注意区分。
求导时注意正负号。
上机实习题三
一、题目:
用Newton法求方程
74+14=0
-28xX
在(0.1,1.9)中的近似根(初始近似值取为区间端点,迭代6次或误差小于0.00001)。
二、解题方法及理论依据:
Newton迭代法是平方收敛于方程f(x)=0在区间[a,b]上的唯一解α,收敛速度较快,循环次数少。
方法应用条件:
ⅰ)f(a)f(b)<0
”ⅱ)f(x)在区间[a,b]上不变号.
'ⅲ)f(x)≠0
'''ⅳ)|f(c)|/b-a≤|f(c)|其中c是a,b中使min[|f(a),f(b)]达到的一个,则对任意时近似值x€0[a,b],Newton迭代过程为
'x=φ(x)=x-f(x)/f(x),k=1,2,3…kk+1kk算法:
令
74?
14,f(0.1)?
0,f(1(fx)?
xx?
28.9)?
06333?
?
16)?
?
7x0(f(x)?
7x112?
xx3252?
?
0)?
?
?
x336x?
42x(x8f)(x?
42?
?
01f(.9)?
f)(1.9?
故以1.9为起点f(x)?
kx?
x?
?
k?
1k?
)(xf?
k?
x?
1.9?
0-8-
数值分析实习报告
.计算程序:
1三、
#includemath.h
main()
{
floatx,y,f,f1;
scanf(%f,&x);
do{
y=x;
*/的表达式定义f(x)f=pow(y,7)-28*pow(y,4)+14;/**/的表达式'(x)f1=7*pow(y,6)-112*pow(y,3);/*定义f*/迭代法x=y-f/f1;/*Newton}
0.00001*/控制误差小于while(fabs(x-y)>=1e-5);/*printf(\
Theresultofthequestionis%f\n,x);
}
.打印结果:
20.1
1.9请输入端点值:
3.030577x=0.845497
四、问题讨论:
轴交点,故也称为切线法,它是平方收敛的,的切线与x是f(x)在点x程序较为简单。
它的几何意义为xkk+1是否为零。
)f′(x此处取x=1.9收敛性较好,要注意判断kk
上机实习题四一、题目:
3241.x?
dxsinx5x3x?
7)(=0.00001)。
用Romberg算法求允许误差(ε1二、解题方法及理论依据:
Romberg)方法求数值积分龙贝格((0)=(b-a)/2*[f(a)+f(b)]T1(l)(l-1)l-1lT=(1/2)*{T+(b-a)/2*∑f[a+(2i-1)*(b-a)/2]}11k-1m(k)(k-1)mT=[4T-T]/(4-1)
mm+1m三、1.计算程序:
#includemath.h
inta=1,b=3;
x2doublef(doublex)/*求f(x)=3x1.4(5x+7)sinx的值*/
-9-
数值分析实习报告
{
doublez;
z=pow(3,x)*pow(x,1.4)*(5*x+7)*sin(x*x);
return(z);
}
ll-1(l)]*/f[a+(2i-1)*(b-a)/2(b-a)/2*∑doubles(intl)/*求T1中的{
externa,b;inti,m;doublez=0.0;
m=pow(2,l-1);
for(i=1;i<=m;i++)z+=f(a+1.0*(2*i-1)/m);
z*=1.0*(b-a)/m;
return(z);
}
main()
{externa,b;
doubleT[20][20];
intm,n,l=0;
T[1][0]=(b-a)/2.0*(f(a)+f(b));
*/)算法龙贝格(Rombergdo/*{l++;
T[1][l]=0.5*(T[1][l-1]+s(l));
n=l-1;
(0)*/求解Tfor(m=2;n>=0;m++,n--)/*lT[m][n]=(pow(4,m-1)*T[m-1][n+1]-T[m-1][n])/(pow(4,m-1)-1.0);
}
while(fabs(T[l][0]-T[l-1][0])>=1e-5);
printf(\
T[%d][0]=%f,l,T[l][0]);
}
.打印结果:
2
T[8][0]=440.536017
四、问题讨论:
(k)外推法,构造新序列,计算新分点的值时,需要复化梯形公式,还要用到RichardsonT此程序较繁,计算1(0)(0)时控制循环。
程序具有广泛的通用性。
|T-T|<εε这些数值个数成倍增加。
应用给出所要求的误差,当l+1l
上机实习题五一、题目:
Runge-Kutta法求解用定步长四阶dy/dt=11dy/dt=y32-10-
数值分析实习报告
-100ydy/dt=1000-1000y332(0)=0y1(0)=0y2(0)=0
y3(0.1),(i=1,2,3)
(0.045),y(0.085),yh=0.0005,打印y(0.025),yiiii二、解题方法及理论依据:
Runge-Kutta解法高阶方程组的Y=Y+(1/6)*(K+2K+2K+K)4n+12n31K=h*F(x,Y)
n1nK=h*F(x+h/2,Y+K/2)1n2nK=h*F(x+h/2,
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 分析 编程 语言 汇总