数值方法作业
- 格式:doc
- 大小:94.00 KB
- 文档页数:3
数值分析上机作业(一)一、算法的设计方案1、幂法求解λ1、λ501幂法主要用于计算矩阵的按模最大的特征值和相应的特征向量,即对于|λ1|≥|λ2|≥.....≥|λn|可以采用幂法直接求出λ1,但在本题中λ1≤λ2≤……≤λ501,我们无法判断按模最大的特征值。
但是由矩阵A的特征值条件可知|λ1|和|λ501|之间必然有一个是最大的,通过对矩阵A使用幂法迭代一定次数后得到满足精度ε=10−12的特征值λ0,然后在对矩阵A做如下的平移:B=A-λ0I由线性代数(A-PI)x=(λ-p)x可得矩阵B的特征值为:λ1-λ0、λ2-λ0…….λ501-λ0。
对B矩阵采用幂法求出B矩阵按模最大的特征值为λ∗=λ501-λ0,所以λ501=λ∗+λ0,比较λ0与λ501的大小,若λ0>λ501则λ1=λ501,λ501=λ0;若λ0<λ501,则令t=λ501,λ1=λ0,λ501=t。
求矩阵M按模最大的特征值λ的具体算法如下:任取非零向量u0∈R nηk−1=u T(k−1)∗u k−1y k−1=u k−1ηk−1u k=Ay k−1βk=y Tk−1u k(k=1,2,3……)当|βk−βk−1||βk|≤ε=10−12时,迭终终止,并且令λ1=βk2、反幂法计算λs和λik由已知条件可知λs是矩阵A 按模最小的特征值,可以应用反幂法直接求解出λs。
使用带偏移量的反幂法求解λik,其中偏移量为μk=λ1+kλ501−λ140(k=1,2,3…39),构造矩阵C=A-μk I,矩阵C的特征值为λik−μk,对矩阵C使用反幂法求得按模最小特征值λ0,则有λik=1λ0+μk。
求解矩阵M按模最小特征值的具体算法如下:任取非零向量u 0∈R n ηk−1= u T (k−1)∗u k−1y k−1=u k−1ηk−1 Au k =y k−1βk =y T k−1u k (k=1,2,3……)在反幂法中每一次迭代都要求解线性方程组Au k =y k−1,当K 足够大时,取λn =1βk 。
问题1:20.给定数据如下表:试求三次样条插值S(x),并满足条件 (1)S`(0.25)=1.0000,S`(0.53)=0.6868; (2)S ’’(0.25)=S ’’(0.53)=0。
分析:本问题是已知五个点,由这五个点求一三次样条插值函数。
边界条件有两种,(1)是已知一阶倒数,(2)是已知自然边界条件。
对于第一种边界(已知边界的一阶倒数值),可写出下面的矩阵方程。
⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡432104321034322110d M M M M M 200020000020022d d d d λμμλμλμλ其中μj =j1-j 1-j h h h +,λi=j1-j j h h h +,dj=6f[x j-1,x j ,x j+1], μn =1,λ0=1对于第一种边界条件d 0=0h 6(f[x 0,x 1]-f 0`),d n =1-n h 6(f`n-f `[x n-1,x n ]) 解:由matlab 计算得:由此得矩阵形式的线性方程组为:⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡ 2.1150-2.4286-3.2667-4.3143-5.5200-M M M M M 25714.00001204286.000004000.026000.0006429.023571.0001243210解得 M 0=-2.0286;M 1=-1.4627;M 2= -1.0333; M 3= -0.8058; M 4=-0.6546S(x)=⎪⎪⎩⎪⎪⎨⎧∈-+-+-∈-+-+-∈-+-+-∈-+-+-]53.0,45.0[x 5.40x 9.1087x 35.03956.8.450-x 1.3637-x .5301.67881- ]45.0,39.0[x 9.30x 11.188x 54.010.418793.0-x 2.2384-x .450(2.87040-]39.0,30.0[x 03.0x 6.9544x 9.30 6.107503.0-x 1.9136-x .3902.708779-]30.0,25.0[x 5.20x 10.9662x 0.3010.01695.20-x 4.8758-x .3006.76209-33333333),()()()(),()()()),()()()(),()()()(Matlab 程序代码如下:function tgsanci(n,s,t) %n代表元素数,s,t代表端点的一阶导。
数值计算方法简明教程第一章1 *1x =1.7; *2x =1.73; *3x =1.732 。
2.3. (1) ≤++)(*3*2*1x x x e r 0.00050; (注意:应该用相对误差的定义去求) (2) ≤)(*3*2*1x x x e r 0.50517; (3) ≤)/(*4*2x x e r 0.50002。
4.设6有n 位有效数字,由6≈2.4494……,知6的第一位有效数字1a =2。
令3)1()1(1*1021102211021)(-----⨯≤⨯⨯=⨯=n n r a x ε 可求得满足上述不等式的最小正整数n =4,即至少取四位有效数字,故满足精度要求可取6≈2.449。
5. 答:(1)*x (0>x )的相对误差约是*x 的相对误差的1/2倍;(2)n x )(* 的相对误差约是*x 的相对误差的n 倍。
6. 根据********************sin 21)(cos 21sin 21)(sin 21sin 21)(sin 21)(c b a c e c b a c b a b e c a c b a a e c b S e r ++≤ =******)()()(tgc c e b b e a a e ++ 注意当20*π<<c 时,0**>>c tgc ,即1*1*)()(--<c tgc 。
则有)()()()(****c e b e a e S e r r r r ++<7.设20=y ,41.1*=y ,δ=⨯≤--2*001021y y 由 δ1*001*111010--≤-=-y y y y ,δ2*111*221010--≤-=-y y y yδ10*991*10101010--≤-=-y y y y即当0y 有初始误差δ时,10y 的绝对误差的绝对值将减小1010-倍。
而11010<<-δ,故计算过程稳定。
丁丽娟《数值计算方法》五章课后实验题答案(源程序都是自己写的,很详细,且保证运行无误)我做的五章数值实验作业题目如下:第二章:1、2、3、4题第三章:1、2题第四章:1、2题第六章:2、3题第八章:1、2题第二章1:(1) 对A进行列主元素三角分解:function [l u]=myfun(A) n=size(A); for k=1:n for i=k:n sum=0; m=k; for j=1:(k-1) sum=sum+A(i,j)*A(j,k); end s(i)=A(i,k)-sum; if abs(s(m))<abs(s(i)) m=i; end end for j=1:n c=A(m,j); A(m,j)=A(k,j); A(k,j)=c; end for j=k:n sum=0; for r=1:(k-1) sum=sum+A(k,r)*A(r,j); end u(k,j)=A(k,j)-sum; A(k,j)=u(k,j); end for i=1:n l(i,i)=1; end for i=(k+1):n sum=0; for r=1:(k-1) sum=sum+A(i,r)*u(r,k); end l(i,k)=(A(i,k)-sum)/u(k,k); A(i,k)=l(i,k); end end 的列主元素三角分解:求A的列主元素三角分解:>>A=[1 1 1 1 1;1 2 3 4 5;1 3 6 10 15;1 4 10 20 35;1 5 15 35 70]; >>[L,U]=myfun(A) 结果:L = 1.0000 0 0 0 0 1.0000 1.0000 0 0 0 1.0000 0.5000 1.0000 0 0 1.0000 0.7500 0.7500 1.0000 0 1.0000 0.2500 0.7500 -1.0000 1.0000 U = 1.0000 1.0000 1.0000 1.0000 1.0000 0 4.0000 14.0000 34.0000 69.0000 0 0 -2.0000 -8.0000 -20.5000 0 0 0 -0.5000 -2.3750 0 0 0 0 -0.2500 (2) 求矩阵的逆矩阵A -1: inv(A) 结果为:ans = 5 -10 10 -5 1 -10 30 -35 19 -4 10 -35 46 -27 6 -5 19 -27 17 -4 1 -4 6 -4 1 (3)检验结果:E=diag([1 1 1 1 1]) A\E ans = 5 -10 10 -5 1 -10 30 -35 19 -4 10 -35 46 -27 6 -5 19 -27 17 -4 1 -4 6 -4 1 2: 程序:程序:function d=myfun(a,b,c,d,n) for i=2:n l(i)=a(i)/b(i-1); a(i)=l(i); u(i)=b(i)-c(i-1)*a(i); b(i)=u(i); y(i)=d(i)-a(i)*d(i-1); d(i)=y(i); end x(n)=d(n)/b(n); d(n)=x(n); for i=(n-1):-1:1 x(i)=(d(i)-c(i)*d(i+1))/b(i); d(i)=x(i); end 求各段电流量程序:求各段电流量程序:for i=2:8 a(i)=-2; end b=[2 5 5 5 5 5 5 5]; c=[-2 -2 -2 -2 -2 -2 -2]; V=220; R=27; d=[V/R 0 0 0 0 0 0 0]; n=8; I=myfun(a,b,c,d,n) 运行程序得:运行程序得:I = 8.1478 4.0737 2.0365 1.0175 0.5073 0.2506 0.1194 0.0477 3:程序:(1)求矩阵A和向量b的matlab程序:function [A b]=myfun(n) for i=1:n X(i)=1+0.1*i; end for i=1:n for j=1:n A(i,j)=X(i)^(j-1); end end for i=1:n b(i)=sum(A(i,:)); end 求n=5时A1,b1及A1的2-条件数程序运行结果如下:条件数程序运行结果如下: n=5;[A1,b1]=myfun(n) A1 = 1.0000 1.1000 1.2100 1.3310 1.4641 1.0000 1.2000 1.4400 1.7280 2.0736 1.0000 1.3000 1.6900 2.1970 2.8561 1.0000 1.4000 1.9600 2.7440 3.8416 1.0000 1.5000 2.2500 3.3750 5.0625 b1 = 6.1051 7.4416 9.0431 10.9456 13.1875 cond2=cond(A1,2)cond2 = 5.3615e+005 条件数程序运行结果如下:求n=10时A2,b2及A2的2-条件数程序运行结果如下:n=10; [A2,b2]=myfun(n) A2 = 1.0000 1.1000 1.2100 1.3310 1.4641 1.6105 1.7716 1.9487 2.1436 2.3579 1.0000 1.2000 1.4400 1.7280 2.0736 2.4883 2.9860 3.5832 4.2998 5.1598 1.0000 1.3000 1.6900 2.1970 2.8561 3.7129 4.8268 6.2749 8.1573 10.6045 1.0000 1.4000 1.9600 2.7440 3.8416 5.3782 7.5295 10.5414 14.7579 20.6610 1.0000 1.5000 2.2500 3.3750 5.0625 7.5938 11.3906 17.0859 25.6289 38.4434 1.0000 1.6000 2.5600 4.0960 6.5536 10.4858 16.7772 26.8435 42.9497 68.7195 1.0000 1.7000 2.8900 4.9130 8.3521 14.1986 24.1376 41.0339 69.7576 118.5879 1.0000 1.8000 3.2400 5.8320 10.4976 18.8957 34.0122 61.2220 110.1996 198.3593 1.0000 1.9000 3.6100 6.8590 13.0321 24.7610 47.0459 89.3872 169.8356 322.6877 1.0000 2.0000 4.0000 8.0000 16.0000 32.0000 64.0000 128.0000 256.0000 512.0000 b2 = 1.0e+003 * 0.0159 0.0260 0.0426 0.0698 0.1133 0.1816 0.2866 0.4451 0.6801 1.0230 cond2=cond(A2,2) cond2 = 8.6823e+011 条件数程序运行结果如下:求n=20时A3,b3及A3的2-条件数程序运行结果如下:n=20; [A3,b3]=myfun(n) A3 = 1.0e+009 * Columns 1 through 10 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Columns 11 through 20 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0001 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0001 0.0001 0.0002 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0001 0.0001 0.0003 0.0005 0.0000 0.0000 0.0000 0.0000 0.0000 0.0001 0.0001 0.0003 0.0006 0.0013 0.0000 0.0000 0.0000 0.0000 0.0001 0.0001 0.0003 0.0007 0.0015 0.0032 0.0000 0.0000 0.0000 0.0001 0.0001 0.0003 0.0006 0.0014 0.0032 0.0075 0.0000 0.0000 0.0000 0.0001 0.0002 0.0005 0.0012 0.0029 0.0070 0.0167 0.0000 0.0000 0.0001 0.0001 0.0004 0.0009 0.0023 0.0058 0.0146 0.0364 0.0000 0.0000 0.0001 0.0002 0.0006 0.0017 0.0044 0.0113 0.0295 0.0766 0.0000 0.0001 0.0002 0.0004 0.0011 0.0030 0.0080 0.0215 0.0581 0.1570 0.0000 0.0001 0.0002 0.0007 0.0018 0.0051 0.0143 0.0400 0.1119 0.3133 0.0000 0.0001 0.0004 0.0010 0.0030 0.0086 0.0250 0.0726 0.2105 0.6103 0.0001 0.0002 0.0005 0.0016 0.0048 0.0143 0.0430 0.1291 0.3874 1.1623 b3 = 1.0e+009 * Columns 1 through 10 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0001 0.0002 0.0004 0.0010Columns 11 through 20 0.0025 0.0059 0.0132 0.0287 0.0606 0.1246 0.2494 0.4874 0.9316 1.7434 cond2=cond(A3,2) cond2 =3.2395e+022 由上述运行结果可知:它们是病态的,而且随着n的增大,矩阵的病态变得严重。
数值计算大作业题目一、非线性方程求根1.题目假设人口随时间和当时人口数目成比例连续增长,在此假设下人口在短期内的增长建立数学模型。
(1)如果令()N t 表示在t 时刻的人口数目,β表示固定的人口出生率,则人口数目满足微分方程()()dN t N t dt β=,此方程的解为0()=tN t N e β; (2)如果允许移民移入且速率为恒定的v ,则微分方程变成()()dN t N t vdt β=+, 此方程的解为0()=+(1)t t vN t N e e βββ-;假设某地区初始有1000000人,在第一年有435000人移入,又假设在第一年年底该地区人口数量1564000人,试通过下面的方程确定人口出生率β,精确到410-;且通过这个数值来预测第二年年末的人口数,假设移民速度v 保持不变。
4350001564000=1000000(1)e e βββ+-2.数学原理采用牛顿迭代法,牛顿迭代法的数学原理是,对于方程0)(=x f ,如果)(x f 是线性函数,则它的求根是很容易的,牛顿迭代法实质上是一种线性化方法,其基本思想是将非线性方程0)(=x f 逐步归结为某种线性方程来求解。
设已知方程0)(=x f 有近似根k x (假定0)(≠'x f ),将函数)(x f 在点k x进行泰勒展开,有.))(()()(⋅⋅⋅+-'+≈k k k x x x f x f x f于是方程0)(=x f 可近似地表示为))(()(=-'+k k x x x f x f这是个线性方程,记其根为1k x +,则1k x +的计算公式为)()(1k k k k x f x f x x '-==+,,,2,1,0⋅⋅⋅=k这就是牛顿迭代法,简称牛顿法。
3.程序设计作出函数的图像,大概估计出根的位置fplot('1000*exp(x)+(435*x)*(exp(x)-1)-1564',[0 3]);grid大概估计出初始值x=0.5function [p1,err,k,y]=newton(f,df,p0,delta,max1) % f 是非线性系数 % df 是f 的微商 % p0是初始值% dalta 是给定允许误差 % max1是迭代的最大次数 % p1是牛顿法求得的方程近似解 % err 是p0误差估计 % k 是迭代次数 p0,feval('f',p0) for k=1:max1p1=p0-feval('f',p0)/feval('df',p0); err=abs(p1-p0); p0=p1;p1,err,k,y=feval('f',p1) if(err<delta)|(y==0), break,endp1,err,k,y=feval('f',p1) endfunction y=f(x)y=1000000*exp(x)+435000*(exp(x)-1)/x-1564000; function y=df(x)y=1000000*exp(x)+435000*(exp(x)/x-(exp(x)-1)/x^2);4.结果分析与讨论newton('f','df',1.2,10^(-4),10) 运行后得出结果 p0 =0.5000p1 =0.1679 err =0.3321 k =1 y =9.2415e+004 p1 =0.1031 err =0.0648 k =2 y =2.7701e+003 p1 =0.1010 err =0.0021 k =3 y =2.6953p1 =0.1010 err =2.0129e-006 k =4 y = 2.5576e-006 ans =0.1010运算后的结果为1010.0=β,通过这个数值来预测第二年年末的人口数,0.10100.1010435000f(t)=1000000(1)0.1010t te e +-t=2时候对于f ()2187945.865x =实践表明,当初始值难以确定时,迭代法就不一定收敛了,因此要根据问题实际背景或者二分法先得一个较好的初始值,然后再进行迭代;再者迭代函数选择不合适的话,采用不动点迭代法也有可能出现不收敛的情况;因此我采用的是牛顿法。
(4) 北京理工大学函大2004-2005学年第1学期计算机科学与技术专业专升本数值计算方法思考题和习题教科书:《科学与工程计算》廖晓钟赖汝编国防工业出版社 2003年版第1 章思考题p26 1,2,3,4,5第1 章习题pp26-27 1,3,4,5,6,11第2 章思考题p66 1,3,6,7,8,9,12.13第2 章习题pp67-68 2,3,4,5,7,11,12,13,14,17,18第3 章思考题p119 1,3,4,5,6,10,18,19第3 章习题pp119-121 1,2,3,4,5,12,13第4 章思考题p144 1,2,3,4,5,7,8第4 章习题pp144-146 1,2,3,4,5,6,7,10,11,12,13第5 章思考题p207 1,2,3,4,5,6,7,9,10,11,12.13第5 章习题pp208-209 1,2,3,4,5,6,7,8,9,10,11,12,13,15第6 章思考题p257 1,2,3,4,5,6,7,8,10,11,12.14第6 章习题pp257-259 1,2,3,4,5,6,7,8,11,12,13,15,16,17,18第7 章思考题p292 1,2,3,4,5,6,8,9第7 章习题pp293-295 1,2,3,4,5,6,7,8,11,12,20作业题第1 章习题pp26-27 1(1),(2),3(3),5,6第2 章习题pp67-68 2,4,5,11,13,17第3 章习题pp119-121 1(1),2(1),5(2),12第4 章习题pp144-146 1(1),2,10,11,12,13第5 章习题pp208-209 1,3,4,7,10,13,,15第6 章习题pp257-259 1(2),3,6(1),12,16第7 章习题pp293-295 1,3,6,11,20数值计算方法复习题第1 章绪论1.说明数值算法的意义,计算机解题步骤和数值算法的特点。
数值分析作业及答案Chap11、写出下列语句的运行结果。
在MA TLAB 上执行它们以验证所得解答。
a=[1 2 3 ;4 5 6 ]’ b=[9;7;5;3;1] c=b(2:4) d=b(4:-1:1) e=sort(b) f=[3,b ’]解:a=635241 b=13579c=357d=9753 e=97531F=[3 9 7 5 3 1] 3、给定一向量:a=[4 -1 2 -8 4 5 -3 -1 6 -7]写一段程序计算a 中正数的和。
运行程序并显示结果。
解:a=[4 -1 2 -8 4 5 -3 -1 6 -7]; s=0;for i=1:length(a) if a(i)>0s=s+a(i); end end s6、编写一个函数M 文件fun_es(x),计算如下函数:230.5sin x y e x x =-其中参数可以为标量,也可以为向量。
在MA TLAB 里键入如下命令检验此函数:fun_es(3) fun_es([1 2 3])解:function y=fun_es(x) y=0.5*exp(x/3)-x.^2.*sin(x);chap21、设0x >,x 的相对误差为δ,求L nx 的误差。
解:Lnx-Lnx*=dLnx=dx/x=δ2、设x 的相对误差为2%,求2x 的相对误差。
解:dLnf(x)=xf ’(x)/f(x)dLnx=4%5、计算球体积要使相对误差限为1%,问度量半径R 时允许的相对误差限是多少?解:dLnf(x)=xf ’(x)/f(x)dLnx=3dLnx=1% dLnx=0.33%9、正方形的边长大约为100cm ,应怎样测量才能使其面积误差不超过1cm 2? 解:s=x 2s-s*=2x(x-x*)=1x-x*=1/(2x)=1/200=0.5*10-2 即测量边的误差不超过0.005cm 10、设212S gt =,假定g 是准确的,而对t 的测量有±0.1秒的误差,证明当t 增加时S 的绝对误差增加,而相对误差却减少。
课程设计课程名称:设计题目:学号:姓名:完成时间:题目一:非线性方程求根 一 摘要非线性方程的解析解通常很难给出,因此非线性方程的数值解就尤为重要。
本实验通过使用常用的求解方法二分法和Newton 法及改进的Newton 法处理几个题目,分析并总结不同方法处理问题的优缺点。
观察迭代次数,收敛速度及初值选取对迭代的影响。
用Newton 法计算下列方程(1) 310x x --= , 初值分别为01x =,00.45x =,00.65x =; (2) 32943892940x x x +-+= 其三个根分别为1,3,98-。
当选择初值02x =时给出结果并分析现象,当6510ε-=⨯,迭代停止。
解:1)采用MATLAB 进行计算;首先定义了Newton 法:function kk=newton(f,df,x0,tol,N)% Newton Method (牛顿法)% The first parameter f is a external function with respect to viable x.(第一个参数也就是本题所用的函数f )% The second parameter df is the first order diffential function of fx.(第二个参数也就是本体所用函数f 的导数方程df ) % x0 is initial iteration point(初值). % tol is the tolerance of the loop (精度).% N is the maximum number of iterations (循环上限). x=x0;f0=eval(f);df0=eval(df); n=0;disp(' [ n xn xn+1 fn+1 ]'); while n<=N x1=x0-f0/df0; x=x1; f1=eval(f); X=[n,x0,x1,f1]; disp(X);if abs(x0-x1)<tolfprintf('The procedure was successful.') kk=X; return else n=n+1; x0=x1;f0=f1;endendif n==N+1fprintf('the method failed after N iterations. '),kk=0;End我们把Newton法存为.m格式的文件;之后我们运行程序:clear;clc;syms xf=x^3-x-1;df=diff(f,x);x=newton(f,df,1,0.0001,50);x会得到一下结果[ n xn xn+1 fn+1 ]0 1.0000 1.5000 0.87501.0000 1.5000 1.0625 -0.86302.0000 1.0625 1.4940 0.8408到第50次迭代时候会出现该问题:47.0000 1.4898 1.0814 -0.816748.0000 1.0814 1.4898 0.816749.0000 1.4898 1.0814 -0.816750.0000 1.0814 1.4898 0.8167the method failed after N iterations.x =0;同样测试x0=0.45、0.65得不出结果,判断出初值离真值太远,所以我们采用牛顿下山法进行计算迭代:我们定义了其中的f函数和df函数,并且分别存为.m格式的文件,其代码如下:f:function y=f(x)y=x^3-x-1;df:function y=df(x)y=3*x^2-1;之后我们定义newton下山法同时也存为.m的程序:function [x,i]=downnewton(f,df,x0,tol)k=0;i=1;disp(' [ n xn xn+1 fn+1 ]'); while(k==0)fx=feval('f',x0);dfx=feval('df',x0);t=0;u=1;while(t==0)dx=-fx/dfx;x1=x0+u*dx;fx1=feval('f',x1);fx0=feval('f',x0);if(abs(fx1)>abs(fx0));u=u/2;elset=1;endendX=[i,x0,x1,fx1];disp(X);if(abs(fx1)<tol)k=1;elsex0=x1;i=i+1;endendx=x1;i=i;end之后带入x0=0.45;downnewton('f','df',0.45,10^(-6))[ n xn xn+1 fn+1 ]1.0000 0.4500 -0.4155 -0.65622.0000 -0.4155 -0.5857 -0.61523.0000 -0.5857 -0.5754 -0.61514.0000 -0.5754 -0.5782 -0.61515.0000 -0.5782 -0.5773 -0.61516.0000 -0.5773 -0.5774 -0.61517.0000 -0.5774 -0.5773 -0.61518.0000 -0.5773 -0.5774 -0.61519.0000 -0.5774 -0.5774 -0.615110.0000 -0.5774 -0.5774 -0.615111.0000 -0.5774 1.3131 -0.049012.0000 1.3131 1.3248 0.000513.0000 1.3248 1.3247 0.0000ans =1.3247带入x0=0.6;downnewton('f','df',0.6,10^(-6))[ n xn xn+1 fn+1 ]1.0000 0.6000 1.1406 -0.65662.0000 1.1406 1.3668 0.18663.0000 1.3668 1.3263 0.00674.0000 1.3263 1.3247 0.00005.0000 1.3247 1.3247 0.0000ans =1.3247带入x0=1;downnewton('f','df',1,10^(-6))[ n xn xn+1 fn+1 ]1.0000 1.0000 1.5000 0.87502.0000 1.5000 1.3478 0.10073.0000 1.3478 1.3252 0.00214.0000 1.3252 1.3247 0.0000ans =1.32472)同样采用Newton下山法:重新定义f、df:f:function y=f(x)y=x^3+94*x^2-389*x+294;df:function y=df(x)y=3*x^2+188*x-389;再带入初值x0=2;downnewton('f','df',2,5*10^(-6))[ n xn xn+1 fn+1 ]1 2 -98 0ans =-98得出x=-98;分析:先画出该函数的图像;x=(-100:.1:100);ezplot('x^3+94*x^2-389*x+294',[-100 100]) 得出该图像如图:-100-80-60-40-20020406080024681012141618x 105xx 3+94 x 2-389 x+294根据牛顿法的几何解释,在x0=2的点做切线,与y 相交,交点的横坐标值为x=-98则结束了该现象。
姓名:潘成功学号:2011412084 学院:机械科学与工程学院专业:机械设计及理论第一次上机作业:三次样条插值
一、题目:
已知函数f(x)=5/(1+x^2),定义区间[-5,5],将定义区间分成10等分。
利用M方法编制三次样条插值函数S(x),满足第一类边界条件:S'(-5)=0.0740,S'(5)= -0.0740,并输出S(2.5)的值。
二、计算
1、编辑三次样条函数cubicsplineinterpolation.m;
2、输入
x=linspace(-5,5,10);f0=0.0740;fn=-0.0740;f=5./(1+x.^2);x0=2.5;
[a,b,c,d,y0]=cubicsplineinterpolation(x,f,f0,fn,x0)
三、结果
第五章 解线性代数方程组的直接法
第二次上机作业:解线性方程组
一、题目:对习题15中的方程组Ax=b ,取初值x (0)=(1,1,1,1,1,1)T ,要求2
)()1(k k x x -+≤10-5 ⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡-=⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣
⎡⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡--------------62605041010014101001400
1100410010141001014654321x x x x x x 二、编写程序,流程图如下:
三、计算结果 k=26 ⎥⎥⎥⎥⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡= 1.9999981.0000001.9999981.0000001.9999981.000000)26(x
姓名:潘成功学号:2011412084 学院:机械科学与工程学院专业:机械设计及理论
第三次上机作业:解非线性方程组
一、题目:解方程组:
(15x1+10x2)/[(40-30x1-10x2)^2×(15-15x1)]=5e-4;
(15x1+10x2)/[(40-30x1-10x2)×(10-10x2)]=4e-2.
二、计算
1、编辑解非线性方程组程序mulNewton.m;
2、输入:
funcMat=[sym('(15*x1+10*x2)-((40-30*x1-10*x2)^2*(15-15*x1))*5e-4')
sym('(15*x1+10*x2)-((40-30*x1-10*x2)*(10-10*x2))*4e-2')];
var=[sym('x1') sym('x2')];
[r,n]=mulNewton([0.2 0.6], funcMat, var, 0.001)
三、结果。