數(shù)值計(jì)算方法程序設(shè)計(jì)

用matlab編寫拉格朗日插值算法的程序并且以(x=-2.00,f(x)=17.00x=0.00,f(x)=1.00x=1.00,f(x)=2.00x=2.00,f(x)=17.00)為數(shù)據(jù)基礎(chǔ)，在整個(gè)插值區(qū)間上采用拉格朗日插值算法計(jì)算f(x=0.6)，寫出程序源代碼，輸出計(jì)算結(jié)果x0=-2.00；x1=0.00；x2=1.00；x3=2.00；y0=17.00;y1=1.00;y2=2.00;y3=17.00;x=0.6y=(x-x1).*(x-x2).*(x-x3)/((x0-x1).*(x0-x2).*(x0-x3))*y0+(x-x0).*(x-x2).*(x-x3)/((x1-x0).*(x1-x2).*(x1-x3))*y1+(x-x0).*(x-x1).*(x-x3)/((x2-x0).*(x2-x1).*(x2-x3))*y2+(x-x0).*(x-x1).*(x-x2)/((x3-x0).*(x3-x1).*(x3-x2))*y3;disp('y=');disp(y);結(jié)果為：x=0.6000y=0.2560追趕法functionx=zhuiganfa%首先說明：追趕法是適用于三對(duì)角矩陣的線性方程組求解的方法，并不適用于其他類型矩陣。%定義三對(duì)角矩陣A的各組成單元。方程為Ax=d%b為A的對(duì)角線元素(1~n)，a為-1對(duì)角線元素(2~n)，c為+1對(duì)角線元素(1~n-1)。TOC\o"1-5"\h\z% A=[2-10 0% -13-2 0% 0-24 -3% 00 -3 5]a=[0-1-2-3];c=[-1-2-3];b=[2345];d=[61-21];n=length(b);u0=0;y0=0;a(1)=0;%“追”的過程L(1)=b(1)-a(1)*u0;y(1)=(d(1)-y0*a(1))/L(1);u(1)=c(1)/L(1);for『2:(n-1)L(i)=b(i)-a(i)*u(i-1);y(i)=(d(i)-y(i-1)*a(i))/L(i);u(i)=c(i)/L(i);endL(n)=b(n)-a(n)*u(n-1);y(n)=(d(n)-y(n-1)*a(n))/L(n);%“趕”的過程x(n)=y(n);fori=(n-1):-1:1x(i)=y(i)-u(i)*x(i+1);end特征向量的計(jì)算，幕法5.2.2冪法的MATLAB程序用冪法計(jì)算矩陣A的主特征值和對(duì)應(yīng)的特征向量的MATLAB主程序function[k,lambda,Vk,Wc]=mifa(A,V0,jd,max1)lambda=0;k=1;Wc=1;,jd=jd*0.1;state=1;V=V0;while((k<=max1)&(state==1))Vk=A*V;[mj]=max(abs(Vk));mk=m;tzw=abs(lambda-mk);Vk=(1/mk)*Vk;Txw=norm(V-Vk);Wc=max(Txw,tzw);V=Vk;lambda=mk;state=0;if(Wc>jd)state=1;endk=k+1;Wc=Wc;endif(Wc<=jd)disp(-請(qǐng)注意：迭代次數(shù)k,主特征值的近似值lambda,主特征向量的近似向量Vk,相鄰兩次迭代的誤差Wc如下：')elsedisp(-請(qǐng)注意：迭代次數(shù)k已經(jīng)達(dá)到最大迭代次數(shù)max1，主特征值的迭代值lambda,主特征向量的迭代向量Vk,相鄰兩次迭代的誤差Wc如下：-)endVk=V;k=k-1;Wc;例5.2.2用幕法計(jì)算下列矩陣的主特征值和對(duì)應(yīng)的特征向量的近似向量，精度e=10-5.并把(1)和(2)輸出的結(jié)果與例5.1.1中的結(jié)果進(jìn)行比較.r123)r122.'-41401B=213C=1-11D=-5130⑵0367;(3)i4-1217；(4)l-102匕7解(1)輸入MATLAB程序>>A=[1-1;24];V0=[1,1]';[k,lambda,Vk,Wc]=mifa(A,V0,0.00001,100),[V,D]=eig(A),Dzd=max(diag(D)),wuD=abs(Dzd-lambda),wuV=V(:,2)./Vk,運(yùn)行后屏幕顯示結(jié)果請(qǐng)注意：迭代次數(shù)k,主特征值的近似值lambda,主特征向量的近似向量Vk,相鄰兩次迭代的誤差Wc如下：k= lambda= Wc=33 3.00000173836804 8.691862856124999e-007Vk= V= wuV=-0.49999942054432 -0.70710678118655 0.44721359549996-0.894428227562941.00000000000000 0.70710678118655 -0.89442719099992-0.89442719099992Dzd=wuD=3 1.738368038406435e-006由輸出結(jié)果可看出，迭代33次，相鄰兩次迭代的誤差叱-8.6919e-007，矩陣A的主特征值的近似值lambda"3.00000和對(duì)應(yīng)的特征向量的近似向量Vk"(-0.50000，1.00000)T,lambda與例5.1.1中A的最大特征值*2=3近似相等，絕對(duì)誤差約為1.73837e-006，

1](,-1)TVk與特征向量XT=匕2 (k2主0)的第1個(gè)分量的絕對(duì)誤差約等于0,第2個(gè)分量的絕對(duì)值相同.由wuV可以看出，氣的特征向量V(:,2)與*的對(duì)應(yīng)分量的比值近似相等因此，用程序mifa.m計(jì)算的結(jié)果達(dá)到預(yù)先給定的精度e=10-5.(2)輸入MATLAB程序>>B=[123;213;336];V0=[1,1,1]';(2)輸入MATLAB程序>>B=[123;213;336];V0=[1,1,1]';[k,lambda,Vk,Wc]=mifa(B,V0,0.00001,100),[V,D]=eig(B),Dzd=max(diag(D)),wuD=abs(Dzd-lambda),wuV=V(:,3)./Vk,運(yùn)行后屏幕顯示結(jié)果請(qǐng)注意：迭代次數(shù)k,主特征值的近似值lambda,主特征向量的近似向量Vk,相鄰兩次迭代的誤差Wc如下：k=3Vk=0.500000000000000.500000000000001.00000000000000V=0.70710678118655-0.70710678118655lambdaWc=0wuV=0.816496580927730.816496580927730.81649658092773Dzd=9wuD0.577350269189630.57735026918963-0.577350269189630.408248290463860.408248290463860.81649658092773(3)輸入MATLAB程序>>C=[122;1-11;4-121];V0=[1,1,1]';[k,lambda,Vk,Wc]=mifa(C,V0,0.00001,100),[V,D]=eig(C),Dzd=max(diag(D)),wuD=abs(Dzd-lambda),Vzd=V(:,1),wuV=V(:,1)./Vk,運(yùn)行后屏幕顯示W(wǎng)c=2.37758124193119wuD=0.90909090909091wuV=0.904534033733350.30151134457778-0.30151134457776請(qǐng)注意：迭代次數(shù)k已經(jīng)達(dá)到最大迭代次數(shù)max1,主特征值的迭代值lambda,主特征向量的迭代向量Vk,相鄰兩次迭代的誤差Wc=2.37758124193119wuD=0.90909090909091wuV=0.904534033733350.30151134457778-0.30151134457776100 0.09090909090910Vzd=0.904534033733290.30151134457776-0.30151134457776DzdVzd=0.904534033733290.30151134457776-0.30151134457776由輸出結(jié)果可見，迭代次數(shù)k已經(jīng)達(dá)到最大迭代次數(shù)max1=100,并且lambda的相鄰兩次迭代的誤差Wc-2.37758>2，由wuV可以看出，lambda的特征向量Vk與真值Dzd的特征向量Vzd對(duì)應(yīng)分量的比值相差較大，所以迭代序列發(fā)散.實(shí)際上，實(shí)數(shù)矩陣C的特征值的近似值為入1T.00000000000001，氣=-i，%=L并且對(duì)應(yīng)的特征向量的近似向量分別為X：=k1(0.90453403373329，0.30151134457776，-0.30151134457776)T，XT=k2(-0.72547625011001,-0.21764287503300-0.07254762501100i,0.58038100008801-0.29019050004400i)T，XT=k3(-0.72547625011001,-0.21764287503300+0.07254762501100i,0.58038100008801+0.29019050004400i)T(k1'0,k2'0,k3'0是常數(shù)).(4)輸入MATLAB程序>>D=[-4140;-5130;-102];V0=[1,1,1]';[k,lambda,Vk,Wc]=mifa(D,V0,0.00001,100),[V,Dt]=eig(D),

Dtzd=max(diag(Dt)),wuDt=abs(Dtzd-lambda),Vzd=V(:,2),wuV=V(:,2)./Vk,運(yùn)行后屏幕顯示結(jié)果請(qǐng)注意：迭代次數(shù)k,主特征值的近似值lambda,主特征向量的近似向量Vk,相鄰兩次迭代的誤差Wc如下：k=lambda= Wc=19 6.00000653949528 6.539523793591684e-006Dtzd= wuDt=6.00000000000000 6.539495284840768e-006Vk=0.797400480535640.71428594783886Vk=0.797400480535640.71428594783886-0.24999918247180Vzd=0.797400480535640.56957177181117-0.19935012013391wuV=0.797400480535640.797400219806180.79740308813370(一)原點(diǎn)位移反冪法的MATLAB主程序1用原點(diǎn)位移反冪法計(jì)算矩陣A的特征值和對(duì)應(yīng)的特征向量的MATLAB主程序1function[k,lambdan,Vk,Wc]=ydwyfmf(A,V0,jlamb,jd,max1)[n,n]=size(A);A1=A-jlamb*eye(n);jd=jd*0.1;RA1=det(A1);ifRA1==0disp('請(qǐng)注意：因?yàn)锳-aE的n階行列式hl等于零，所以A-aE不能進(jìn)行LU分解.')returnendlambda=0;ifRA1~=0forp=1:nh(p)=det(A1(1:p,1:p));endhl=h(1:n);fori=1:nifh(1,i)==0disp(-請(qǐng)注意：因?yàn)锳-aE的r階主子式等于零，所以A-aE不能進(jìn)行LU分解.')returnendendifh(1,i)?=0disp(-請(qǐng)注意：因?yàn)锳-aE的各階主子式都不等于零，所以A-aE能進(jìn)行LU分解.')k=1;Wc=1;state=1;Vk=V0;while((k<=max1)&(state==1))[LU]=lu(A1);Yk=L\Vk;Vk=U\Yk;[mj]=max(abs(Vk));mk=m;Vk1=Vk/mk;Yk1=L\Vk1;Vk1=U\Yk1;[mj]=max(abs(Vk1));mk1=m;Vk2=(1/mk1)*Vk1;tzw1=abs((mk-mk1)/mk1);tzw2=abs(mk1-mk);Txw1=norm(Vk)-norm(Vk1);Txw2=(norm(Vk)-norm(Vk1))/norm(Vk1);Txw=min(Txw1,Txw2);tzw=min(tzw1,tzw2);Vk=Vk2;mk=mk1;Wc=max(Txw,tzw);Vk=Vk2;mk=mk1;state=0;if(Wc>jd)state=1;endk=k+1;%Vk=Vk2,mk=mk1,endif(Wc<=jd)disp('A-aE的秩R(A-aE)和各階順序主子式值hl、迭代次數(shù)k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相鄰兩次迭代的誤差Wc如

elsedisp('A-aE的秩R(A-aE)和各階順序主子式值hl、迭代次數(shù)k已經(jīng)達(dá)到最大迭代次數(shù)maxi，按模最小特征值的迭代值lambda,特征向量的迭代向量Vk,相鄰兩次迭代的誤差Wc如下：')endhl,RAiendend[V,D]=eig(A,'nobalance'),Vk;k=k-1;Wc;lambdan=jlamb+1/mk1;… 入.例5.3.2用原點(diǎn)位移反幕法的迭代公式(5.28)，根據(jù)給定的下列矩陣的特征值n的初始值值n，計(jì)算與七對(duì)應(yīng)的特征向量Xn的近似向量,精確到0.0001."1-101r-11215-24-2r1-1)2583?(1)3-12J人=0.2(2),2 ；(2)<24/人=2.001悠)，2 ；(3)、153-3人=8.26，3 .解(1)輸入MATLAB程序>>A=[1-10;-24-2;0-12];V0=[1,1,1]';[k,lambda,Vk,Wc]=ydwyfmf(A,V0,0.2,0.0001,10000)運(yùn)行后屏幕顯示結(jié)果請(qǐng)注意：因?yàn)锳-aE的各階主子式都不等于零，所以A-aE能進(jìn)行LU分解.A-aE的秩A-aE的秩R(A-aE)和各階順序主子式值hl、迭代次數(shù)k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相鄰兩次迭代的誤差Wc如下:30.23841.0213e-0070.80001.04000.2720Vk=V=D=1.0000-0.2424-1.0000 -0.57075.12490 00.76161.0000-0.7616 0.363300.2384 00.4323-0.3200-0.4323 1.000000 1.6367k=lambda=Wc=hl=(2)輸入MATLAB程序>>A=[1-1;24];V0=[20,1]';[k,lambda,Vk,Wc]=ydwyfmf(A,V0,2.001,0.0001,100)運(yùn)行后屏幕顯示結(jié)果請(qǐng)注意：因?yàn)锳-aE的各階主子式都不等于零，所以A-aE能進(jìn)行LU分解.A-aE的秩R(A-aE)和各階順序主子式值hl、迭代次數(shù)k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相鄰兩次迭代的誤差Wc如下：k=lambda=Wc=hl=22.00205.1528e-007-1.0010-0.0010Vk=V=D=1.0000-1.0000 0.500020-1.00001.0000 -1.000003(3)輸入MATLAB程序>>A=[-11215;2583;153-3];V0=[1,1,-1]';[k,lambdan,Vk,Wc]=ydwyfmf(A,V0,8.26,0.0001,100)運(yùn)行后屏幕顯示結(jié)果請(qǐng)注意：因?yàn)锳-aE的各階主子式都不等于零，所以A-aE能進(jìn)行LU分解.A-aE的秩R(A-aE)和各階順序主子式值hl、迭代次數(shù)k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相鄰兩次迭代的誤差Wc如下：k=lambdan=Wc=hl=28.26406.9304e-008-19.2600-961.9924-6.1256Vk=V=D=-0.76920.79280.60810.0416-22.5249 000.09120.0030-0.07210.99740 8.26400-1.0000-0.60950.79060.05900 058.2609

f011一5]A=-217-7例5.3.3用原點(diǎn)位移反幕法的迭代公式(5.28),計(jì)算"―426-10J的分別對(duì)應(yīng)工蛀行估入rX=1.001入~X=2.001"“"—4.°01m蛀什閂旦XXX,,于特征值11 , 2 2 ， 3 3 的特征向里1, 2， 3的近似向量，相鄰迭代誤差為0.001.將計(jì)算結(jié)果與精確特征向量比較.解(1)計(jì)算特征值A(chǔ)rR—1.001對(duì)應(yīng)的特征向量X1的近似向量.輸入MATLAB程序>>A=[011-5;-217-7;-426-10];V0=[1,1,1]';[k,lambda,Vk,Wc]=ydwyfmf(A,V0,1.001,0.001,100),[V,D]=eig(A);Dzd=min(diag(D)),wuD=abs(Dzd-lambda),VD=V(:,1),wuV=V(:,1)./Vk,運(yùn)行后屏幕顯示結(jié)果請(qǐng)注意：因?yàn)锳-aE的各階主子式都不等于零，所以A-aE能進(jìn)行Lu分解.A-aE的秩R(A-aE)和各階順序主子式值hl、迭代次數(shù)k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相鄰兩次迭代的誤差Wc如下：hl=-1.00100000000000 5.98500100000000 -0.00299600100000RA1=-0.00299600100000k= lambdaRA1=-0.002996001000005 1.00200000000000Vk=-0.50000000000000Vk=-0.50000000000000-0.50000000000000-1.00000000000000Wc=1.378794763695562e-009VD=-0.40824829046386-0.40824829046386-0.81649658092773Dzd=1.00000000000000wuV=0.816496580927730.816496580927730.81649658092773wuD=0.00200000000000~從輸出的結(jié)果可見，迭代5次，特征向量X1的近似向量X1的相鄰兩次迭代的誤差~WcR1.379e-009由wuV可以看出，X1=Vk與VD的對(duì)應(yīng)分量的比值相等.特征值X1的近?似值lambdaR1.002與初始值入1=1.001的絕對(duì)誤差為0.001，而與X1的絕對(duì)誤差為0.002,其中X1—(-0.50000000000000…，-0.50000000000000…，1.00000000000000…)T?X=(-0.50000000000000,-0.50000000000000,1.00000000000000)T1.(2)計(jì)算特征值X2rX2=2.001對(duì)應(yīng)特征向量X2的近似向量.輸入MATLAB程序>>A=[011-5;-217-7;-426-10];V0=[1,1,1]';[k,lambda,Vk,Wc]=ydwyfmf(A,V0,2.001,0.001,100),[V,D]=eig(A);WD=lambda-D(2,2),VD=V(:,2),wuV=V(:,2)./Vk,運(yùn)行后屏幕顯示結(jié)果請(qǐng)注意：因?yàn)锳-aE的各階主子式都不等于零，所以A-aE能進(jìn)行Lu分解.A-aE的秩R(A-aE)和各階順序主子式值hl、迭代次數(shù)k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相鄰兩次迭代的誤差Wc如下：hl=-2.00100000000000 -8.01299900000000 0.00200099900000k=Wc= lambda= WD=2 3.131363162302120e-007 2.002000000000160.00200000000016Vk=-0.24999999999999-0.49999999999999-1.00000000000000VDVk=-0.24999999999999-0.49999999999999-1.00000000000000VD=0.218217890235990.436435780471980.87287156094397wuV=-0.87287156094401-0.87287156094398-0.87287156094397~從輸出的結(jié)果可見，迭代2次，特征向量X2的近似向量X2的相鄰兩次迭代的誤差g=3.131e-007，X2與X2的對(duì)應(yīng)分量的比值近似相等.特征值*2的近似值lambda-2.002與初始值*2=2.001的絕對(duì)誤差約為0.001，而lambda與*2的絕對(duì)誤差約為0.002，其中~X=(-0.24999999999999,-0.49999999999999,-1.00000000000000)tX2=(-0.24999999999999…，-0.50000000000000…，-1.00000000000000…)t(3)計(jì)算特征值*3°*3=4，001對(duì)應(yīng)特征向量X3的近似向量.輸入MATLAB程序>>A=[011-5;-217-7;-426-10];V0=[1,1,1]';[k,lambda,Vk,Wc]=ydwyfmf(A,V0,4.001,0.001,100)[V,D]=eig(A);WD=lambda-max(diag(D)),VD=V(:,3),wuV=V(:,3)./Vk,運(yùn)行后屏幕顯示結(jié)果請(qǐng)注意：因?yàn)锳-aE的各階主子式都不等于零，所以A-aE能進(jìn)行Lu分解.A-aE的秩R(A-aE)和各階順序主子式值hl、迭代次數(shù)k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相鄰兩次迭代的誤差Wc如下：hl=-4.00100000000000 -30.00899900000000 -0.00600500099999WD=0.00199999999990k=lambda= WcWD=0.001999999999902 4.00199999999990 1.996084182914842e-007Vk=0.400000000000010.60000000000001Vk=0.400000000000010.600000000000011.00000000000000VD=-0.32444284226153-0.48666426339229-0.81110710565381wuV=-0.81110710565380-0.81110710565381-0.81110710565381~從輸出的結(jié)果可見，迭代2次，特征向量X3的近似向量X3的相鄰兩次迭代的誤差~Wc°1.996e-007，X3與X3的對(duì)應(yīng)分量的比值近似相等.特征值*3的近似值lambda°4-002與初始值*2=4'001的絕對(duì)誤差近似為0.001，而lambda與*3的絕對(duì)誤差約為0.002,其中X=(3 (-0.40000000000000，-0.60000000000000，-1.00000000000000)T,~X3=(0.4000000000001,0.60000000000001,1.00000000000000)T(二)原點(diǎn)位移反冪法的MATLAB主程序2用原點(diǎn)位移反冪法計(jì)算矩陣A的特征值和對(duì)應(yīng)的特征向量的MATLAB主程序2function[k,lambdan,Vk,Wc]=wfmifa1(A,V0,jlamb,jd,max1)[n,n]=size(A);jd=jd*0.1;A1=A-jlamb*eye(n);nA1=inv(A1);lambda1=0;k=1;Wc=1;state=1;U=V0;while((k<=max1)&(state==1))Vk=A1\U;[mj]=max(abs(Vk));mk=m;Vk=(1/mk)*Vk;Vk1=A1\Vk;[m1j]=max(abs(Vk1));mk1=m1,Vk1=(1/mk1)*Vk1;U=Vk1,Txw=(norm(Vk1)-norm(Vk))/norm(Vk1);tzw=abs((lambda1-mk1)/mk1);Wc=max(Txw,tzw);lambda1=mk1;state=0;if(Wc>jd)state=1;endk=k+1;endif(Wc<=jd)

disp('請(qǐng)注意迭代次數(shù)k,特征值的近似值lambda,對(duì)應(yīng)的特征向量的近似向量Vk,相鄰兩次迭代的誤差Wc如下：')elsedisp(-請(qǐng)注意迭代次數(shù)k已經(jīng)達(dá)到最大迭代次數(shù)maxi,特征值的近似值lambda,對(duì)應(yīng)的特征向量的近似向量Vk,相鄰兩次迭代的誤差Wc如下：-)end[V,D]=eig(A,'nobalance'),Vk=U;k=k-1;Wc;lambdan=jlamb+1/mk;例5.3.4用原點(diǎn)位移反幕法的迭代公式(5.27)，計(jì)算例題5.3.3,并且將這兩個(gè)例題的計(jì)算結(jié)果進(jìn)行比較.再用兩種原點(diǎn)位移反幕法的MATLAB主程序，求「0.99999999999997對(duì)應(yīng)的特征向量.解(1)計(jì)算特征值X1^^1=1.001對(duì)應(yīng)特征向量X1的近似向量.輸入MATLAB程序>>A=[011-5;-217-7;-426-10];V0=[1,1,1]';[k,lambda,Vk,Wc]=wfmifa1(A,V0,1.001,0.001,100)運(yùn)行后屏幕顯示結(jié)果請(qǐng)注意迭代次數(shù)k,特征值的近似值lambda,對(duì)應(yīng)的特征向量的近似向量Vk,相鄰兩次迭代的誤差Wc如下：k=lambda= Wc=5 1.00200000000138 1.376344154436924e-006Vk’= -0.50000000000000 -0.50000000000000-1.00000000000000同理可得，另外與兩個(gè)特征值對(duì)應(yīng)的特征向量.⑵再用兩種原點(diǎn)位移反冪法的MATLAB主程序，求%=0.99999999999997對(duì)應(yīng)的特征向量.輸入MATLAB程序>>A=[011-5;-217-7;-426-10];V0=[1,1,1]';[k,lambda,Vk,Wc]=ydwyfmf(A,V0,0.99999999999997,0.001,100)運(yùn)行后屏幕顯示結(jié)果請(qǐng)注意：因?yàn)锳-aE的各階主子式都不等于零，所以A-aE能進(jìn)行LU分解.A-aE的秩R(A-aE)和各階順序主子式值hl、迭代次數(shù)k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相鄰兩次迭代的誤差Wc如下：hl=-0.99999999999997 6.00000000000045 0.00000000000010Vk=Vk=0.500000000000000.500000000000001.00000000000000Wc=4.317692037236759e-013RA1=1.039168751049192e-013

k=2

lambda=1.00000000000000輸入MATLAB程序>>A=[011-5;-217-7;-426-10];V0=[1,1,1]';[k,lambda,Vk,Wc]=wfmifa1(A,V0,0.99999999999997,0.001,100)運(yùn)行后屏幕顯示結(jié)果請(qǐng)注意迭代次數(shù)k,特征值的近似值lambda,對(duì)應(yīng)的特征向量的近似向量Vk,相鄰兩次迭代的誤差Wc如下：Vk= k=0.50000000000000 30.50000000000000 lambda=1.00000000000000 1.00000000000000Wc=5.412337245047640e-0165.4雅可比(Jacobi)方法及其MATLAB程序5.4.3雅可比方法的MATLAB程序用雅可比方法計(jì)算對(duì)稱矩陣A的特征值和對(duì)應(yīng)的特征向量的MATLAB主程序function[k,Bk,V,D,Wc]=jacobite(A,jd,max1)[n,n]=size(A);Vk=eye(n);Bk=A;state=1;k=0;P0=eye(n);Aij=abs(Bk-diag(diag(Bk)));[m1i]=max(Aij);[m2j]=max(m1);i=i(j);while((k<=max1)&(state==1))k=k+1,aij=abs(Bk-diag(diag(Bk)));[m1i]=max(abs(aij));[m2j]=max(m1);i=i(j),j,Aij=(Bk-diag(diag(Bk)));mk=m2*sign(Aij(i,j)),Wc=m2,Dk=diag(diag(Bk));Pk=P0;c=(Bk(j,j)-Bk(i,i))/(2*Bk(i,j)),t=sign(c)/(abs(c)+sqrt(1+c八2)),pii=1/(sqrt(1+t八2)),pij=t/(sqrt(1+t八2)),Pk(i,i)=pii;Pk(i,j)=pij;Pk(j,j)=pii;Pk(j,i)=-pij;Pk,B1=Pk'*Bk;B2=B1*Pk;Vk=Vk*Pk,Bk=B2,if(Wc>jd)state=1;elsereturnendPk;Vk;Bk=B2;Wc;endif(k>max1)disp(-請(qǐng)注意迭代次數(shù)k已經(jīng)達(dá)到最大迭代次數(shù)max1，迭代次數(shù)k,對(duì)稱矩陣Bk,以特征向量為列向量的矩陣V,特征值為對(duì)角元的對(duì)角矩陣D如下：')elsedisp(-請(qǐng)注意迭代次數(shù)k,對(duì)稱矩陣Bk,以特征向量為列向量的矩陣V,特征值為對(duì)角元的對(duì)角矩陣D如下：')endWc;k=k;V=Vk;Bk=B2;D=diag(diag(Bk));[V1,D1]=eig(A,'nobalance')5.常微分方程數(shù)值解法用matlab解微分方程組dx/dt=x-y-x(xA2+yA2)dy/dt=x+y-y(xA2+yA2)x(0)=2y(0)1解析解：[x,y]=dsolve('Dx=x-y-x*(xA2+yA2)','Dy=x+y-y*(xA2+yA2)','x(0)=2','y(0)=1')得到的結(jié)果是解析解沒有找到。用數(shù)值解。在Matlab下輸入：edit，然后將下面兩行百分號(hào)之間的內(nèi)容，復(fù)制進(jìn)去，保存%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%functiony=zhidao_rk4_5(t,x)%x,y變量分別用x(1),x(2)表示y=[x(1)-x(2)-x(1)*(x(1)A2+x(2)A2);x(1)+x(2)-x(2)*(x(1)A2+x(2)A2)];%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%在Matlab下面輸入：t_end=10;x0=[2;1];[t,x]=ode45('zhidao_rk4_5',[0,t_end],x0);plot(t,x);legend('x','y');xlabel('t');figure;plot(x(:,1),x(:,2));xlabel('x');ylabel('y');6.復(fù)化梯形公式用復(fù)化梯形公式求解sinx積分，積分區(qū)間0—pi建立Trapezoid.m文件function[I,step]=Trapezoid(f,a,b,eps)%f為函數(shù)，a為積分上限，b為積分下限，eps為積分精度，step為劃分區(qū)間個(gè)數(shù)if(nargin==3)eps=1.0e-4;endn=1;h=(b-a)/2;I1=2;I2=(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym(f)),b))/h;whileabs(I2-I1)>epsn=n+1;h=(b-a)/n;I1=I2;I2=0;fori=0:n-1 %第n次的復(fù)化梯形公式積分x=a+h*i;%i=0和n-1時(shí)，分別代表積分區(qū)間的左右端點(diǎn)x1=x+h;I2=I2+(h/2)*(subs(sym(f),findsym(sym(f)),x)+...subs(sym(f),findsym(sym(f)),x1)); %公式endendI=I2;step=n;輸入：[q,s]=Trapezoid('sin(x)',0,pi) %程序默認(rèn)精度為1e-4輸出：q=1.9985;s=33若輸入[q,s]=Trapezoid(‘sin(x)’，0,6pi,1e輸出：q=1.9999;s=157.復(fù)化梯形公式求積分(matlab)2007-11-2719:46程序function[I,n,Ichain]=computT(a,b,errorBound,dNum)%復(fù)化梯形公式求積分%調(diào)用格式：[I,n,Ichain]=computT(a,b,errorBound,dNum)%輸入4：%a:積分下限b：積分上限errorBound:輸出結(jié)果的精度dNum:區(qū)間初始分點(diǎn)數(shù)%輸出3：% I:積分近似值n：最終區(qū)間分點(diǎn)數(shù)Ichain-迭代過程所有值%被積函數(shù)做成函數(shù)文件f(x)n=dNum/2;hn=(b-a)/n;k=1;Tn=hn*(f(a)/2+sumf(a,hn,n-1)+f(b)/2);h2n=hn/2;T2n=computT2n(a,Tn,h2n,n);error=abs(Tn-T2n);whileerror>errorBoundh2n=h2n/2;n=n*2;T2ns=computT2n(a,T2n,h2n,n);Ichain(k)=T2ns;k=k+1;error=abs(T2ns-T2n);T2n=T2ns;endI=T2n;% functionout=computT2n(a,Tn,h2n,n)out=Tn/2+h2n*sumf2(a,h2n,n);out=0;fori=1:n1out=out+f(a+i*hn);endfunctionout=sumf2(a,hn,n1)out=0;fori=1:n1out=out+f(a+(2*i-1)*hn);en