>> norm(X)/norm(x1) ans =
6.3276
6.3276 即为?x >> norm(a) ans =
0.2244 >> norm(A) ans =
30.2887
>> norm(a)/norm(A) %求?Aans =
0.0074
所以?A222x2
A2
A2=0.0074
||A?1||||?A||||x||inv(A) %求A的逆矩阵,下求
1?||A?1(?A)||>> d=inv(A);
>> norm(d)*norm(a)*norm(x) ans =
288.4858
>> 1-norm(d*(a)) ans =
-12.3693
>> 288.4858/ -12.3693 ans =
-23.3227
||A?1||||?A||||x||所以=-23.322
1?||A?1(?A)||>>norm(X) ans =
0.0074
||A?1||||?A||||x||所以?x2>
1?||A?1(?A)||(3)改变?A,取a2=[0 0.002 0.001 0.003;0.001 0.004 0 0.001;0.003 -0.004 -0.001 0;-0.001 -0.002 0 -0.003] B2=a2+A; C2=[B b] b=[32;23;33;31]
>> rank(B2)
ans =
4
>> rank(C2)
ans =
4
rank(B2)= rank(C2)
所以扰动矩阵有唯一解 >> x2=B2\\b x2 =
1.0649 0.8893
>> norm(X)/norm(x1) %求?x2x2
>> norm(X2)/norm(x1) ans =
6.3276
所以?x2x2=6.3276
>> norm(a2)
ans =
0.0071
>> norm(a)/norm(A) %求?A2A2
>> norm(a2)/norm(A)
ans =
2.3336e-004
所以?A2A2=2.3336e-004
1.0272
0.9859
>> x=B\\b;x1=[1;1;1;1];X=x-x1 %求?x(设X= ?x) X2 =
-6.4761 10.4934 -2.4292 1.4838
norm(X2) %求 ?x2 >> norm(X2)
ans =
12.6552
12.6552就是?x2
norm(d)*norm(a2)*norm(x2)
ans =
9.0875
> 1-norm(d*(a2))
ans =
0.6943 norm(X2)
ans =
12.6552
9.0875/0.6943
ans =
13.0887
||A?1||||?A||||x||所以= 13.0887 ?11?||A(?A)||||A?1||||?A||||x||?x2< ?11?||A(?A)||
三、解三对角线性方程组的追赶法及其应用
⑴编写解三对角线性方程组的追赶法的通用程序,并应用于方程组
?2?1000??x1??1???12?100??x??0?T???2????52111??0?12?10??x3???0?,检验程序的正确性;(解为x??,,,,?)
?63236????????00?12?1??x4??0???000?12????0???x5????d2udu?u?ex?3sinx,0?x??1?2?⑵求微分方程边值问题?dx的数值解(取步长h?),dx128?u(0)??2,u(?)?e??3?并与精确解比较(精确解为u?1)。 1?xu?2ui?ui?1duui?1?ui?1d2u,|?,i?1,2,?,n?1 说明:离散化微分方程时,2|x?xi?i?1x?xdxh2dxi2h
解:clear all;
a=[2,2,2,2,2]; b=[-1,-1,-1,-1]; c=[-1,-1,-1,-1]; r=[1,0,0,0,0]; n=length(a); b=[0,b]; u(1)=r(1)/a(1); v(1)=c(1)/a(1); for k=2:n-1
u(k)=(r(k)-u(k-1)*b(k))/(a(k)-v(k-1)*b(k)); v(k)=c(k)/(a(k)-v(k-1)*b(k)); end
u(n)=(r(n)-u(n-1)*b(n))/(a(n)-v(n-1)*b(n)); x(n)=u(n); for k=n-1:-1:1
x(k)=u(k)-v(k)*x(k+1); end
fprintf('èy????·?3ì×éμ??a?a\\n')
for k=1:n fprintf('x()=.8f\\n',k,x(k)) end
>> zhuiganfa %调用追赶法 三对角方程组的解为 x(1)=0.83333333 x(2)=0.66666667 x(3)=0.50000000 x(4)=0.33333333 x(5)=0.16666667
?52111?和结果x??,,,,?很好地吻合。
?63236?
四、公元1225年,比萨的数学家Leonardo研究了方程x?2x?10x?20?0,得到一个根x*?1.368808107,没有人知道他用什么方法得到这个值。对于这个方程,分别用下列
32T20?2xk2?xk320方法:⑴迭代法xk?1?2;⑵迭代法xk?1?;⑶对⑴的Steffensen
xk?2xk?1010加速方法;⑷对⑵的Steffensen加速方法;⑸Newton法。求方程的根(可取x0?1),计算到Leonardo所得到的准确度。
解:由题意编写m文件如下
function[x0,k,er,x]=diedai(g,x0,wucha,max) %g是给定的迭代函数
%x0是给定的初始值 %wucha是规定的误差范围 %max是所应许的最大迭代次数 %k是迭代次数加1 %x是不动点近似值 %x(x1,x2...,xn) X(1)=x0; for k=2:max
X(k)=feval('g',X(k-1)); k,er=abs(X(k)-X(k-1)) x=X(k); if(er disp('超出迭代次数'); end end 其中定义的g函数是 function y=g(x); y=20/(x^2+2*x+10); 在命令窗口中输入>> diedai('g',1,10^(-9),15) k = 15 er = 1.410125245193683e-005 超出迭代次数 X = Columns 1 through 3 1.00000000000000 1.53846153846154 1.29501915708812 Columns 4 through 6 1.40182530944860 1.35420939040429 1.37529809248738 Columns 7 through 9 1.36592978817065 1.37008600340182 1.36824102361284 Columns 10 through 12 1.36905981200748 1.36869639755552 1.36885768862873