数值计算方法程序设计

用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)为数据基础，在整个插值区间上采用拉格朗日插值算法计算f(x=0.6)，写出程序源代码，输出计算结果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);结果为：x=0.6000y=0.2560追赶法functionx=zhuiganfa%首先说明：追赶法是适用于三对角矩阵的线性方程组求解的方法，并不适用于其他类型矩阵。%定义三对角矩阵A的各组成单元。方程为Ax=d%b为A的对角线元素(1~n)，a为-1对角线元素(2~n)，c为+1对角线元素(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特征向量的计算，幕法5.2.2幂法的MATLAB程序用幂法计算矩阵A的主特征值和对应的特征向量的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(-请注意：迭代次数k,主特征值的近似值lambda,主特征向量的近似向量Vk,相邻两次迭代的误差Wc如下：')elsedisp(-请注意：迭代次数k已经达到最大迭代次数max1，主特征值的迭代值lambda,主特征向量的迭代向量Vk,相邻两次迭代的误差Wc如下：-)endVk=V;k=k-1;Wc;例5.2.2用幕法计算下列矩阵的主特征值和对应的特征向量的近似向量，精度e=10-5.并把(1)和(2)输出的结果与例5.1.1中的结果进行比较.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,运行后屏幕显示结果请注意：迭代次数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由输出结果可看出，迭代33次，相邻两次迭代的误差叱-8.6919e-007，矩阵A的主特征值的近似值lambda"3.00000和对应的特征向量的近似向量Vk"(-0.50000，1.00000)T,lambda与例5.1.1中A的最大特征值*2=3近似相等，绝对误差约为1.73837e-006，

1](,-1)TVk与特征向量XT=匕2 (k2主0)的第1个分量的绝对误差约等于0,第2个分量的绝对值相同.由wuV可以看出，气的特征向量V(:,2)与*的对应分量的比值近似相等因此，用程序mifa.m计算的结果达到预先给定的精度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,运行后屏幕显示结果请注意：迭代次数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,运行后屏幕显示Wc=2.37758124193119wuD=0.90909090909091wuV=0.904534033733350.30151134457778-0.30151134457776请注意：迭代次数k已经达到最大迭代次数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由输出结果可见，迭代次数k已经达到最大迭代次数max1=100,并且lambda的相邻两次迭代的误差Wc-2.37758>2，由wuV可以看出，lambda的特征向量Vk与真值Dzd的特征向量Vzd对应分量的比值相差较大，所以迭代序列发散.实际上，实数矩阵C的特征值的近似值为入1T.00000000000001，气=-i，%=L并且对应的特征向量的近似向量分别为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是常数).(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,运行后屏幕显示结果请注意：迭代次数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(一)原点位移反幂法的MATLAB主程序1用原点位移反幂法计算矩阵A的特征值和对应的特征向量的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('请注意：因为A-aE的n阶行列式hl等于零，所以A-aE不能进行LU分解.')returnendlambda=0;ifRA1~=0forp=1:nh(p)=det(A1(1:p,1:p));endhl=h(1:n);fori=1:nifh(1,i)==0disp(-请注意：因为A-aE的r阶主子式等于零，所以A-aE不能进行LU分解.')returnendendifh(1,i)〜=0disp(-请注意：因为A-aE的各阶主子式都不等于零，所以A-aE能进行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、迭代次数k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相邻两次迭代的误差Wc如

elsedisp('A-aE的秩R(A-aE)和各阶顺序主子式值hl、迭代次数k已经达到最大迭代次数maxi，按模最小特征值的迭代值lambda,特征向量的迭代向量Vk,相邻两次迭代的误差Wc如下：')endhl,RAiendend[V,D]=eig(A,'nobalance'),Vk;k=k-1;Wc;lambdan=jlamb+1/mk1;… 入.例5.3.2用原点位移反幕法的迭代公式(5.28)，根据给定的下列矩阵的特征值n的初始值值n，计算与七对应的特征向量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)运行后屏幕显示结果请注意：因为A-aE的各阶主子式都不等于零，所以A-aE能进行LU分解.A-aE的秩A-aE的秩R(A-aE)和各阶顺序主子式值hl、迭代次数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)运行后屏幕显示结果请注意：因为A-aE的各阶主子式都不等于零，所以A-aE能进行LU分解.A-aE的秩R(A-aE)和各阶顺序主子式值hl、迭代次数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)运行后屏幕显示结果请注意：因为A-aE的各阶主子式都不等于零，所以A-aE能进行LU分解.A-aE的秩R(A-aE)和各阶顺序主子式值hl、迭代次数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用原点位移反幕法的迭代公式(5.28),计算"―426-10J的分别对应工蛀行估入rX=1.001入~X=2.001"“"—4.°01m蛀什闩旦XXX,,于特征值11 , 2 2 ， 3 3 的特征向里1, 2， 3的近似向量，相邻迭代误差为0.001.将计算结果与精确特征向量比较.解(1)计算特征值ArR—1.001对应的特征向量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,运行后屏幕显示结果请注意：因为A-aE的各阶主子式都不等于零，所以A-aE能进行Lu分解.A-aE的秩R(A-aE)和各阶顺序主子式值hl、迭代次数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~从输出的结果可见，迭代5次，特征向量X1的近似向量X1的相邻两次迭代的误差~WcR1.379e-009由wuV可以看出，X1=Vk与VD的对应分量的比值相等.特征值X1的近〜似值lambdaR1.002与初始值入1=1.001的绝对误差为0.001，而与X1的绝对误差为0.002,其中X1—(-0.50000000000000…，-0.50000000000000…，1.00000000000000…)T〜X=(-0.50000000000000,-0.50000000000000,1.00000000000000)T1.(2)计算特征值X2rX2=2.001对应特征向量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,运行后屏幕显示结果请注意：因为A-aE的各阶主子式都不等于零，所以A-aE能进行Lu分解.A-aE的秩R(A-aE)和各阶顺序主子式值hl、迭代次数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~从输出的结果可见，迭代2次，特征向量X2的近似向量X2的相邻两次迭代的误差g=3.131e-007，X2与X2的对应分量的比值近似相等.特征值*2的近似值lambda-2.002与初始值*2=2.001的绝对误差约为0.001，而lambda与*2的绝对误差约为0.002，其中~X=(-0.24999999999999,-0.49999999999999,-1.00000000000000)tX2=(-0.24999999999999…，-0.50000000000000…，-1.00000000000000…)t(3)计算特征值*3°*3=4，001对应特征向量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,运行后屏幕显示结果请注意：因为A-aE的各阶主子式都不等于零，所以A-aE能进行Lu分解.A-aE的秩R(A-aE)和各阶顺序主子式值hl、迭代次数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~从输出的结果可见，迭代2次，特征向量X3的近似向量X3的相邻两次迭代的误差~Wc°1.996e-007，X3与X3的对应分量的比值近似相等.特征值*3的近似值lambda°4-002与初始值*2=4'001的绝对误差近似为0.001，而lambda与*3的绝对误差约为0.002,其中X=(3 (-0.40000000000000，-0.60000000000000，-1.00000000000000)T,~X3=(0.4000000000001,0.60000000000001,1.00000000000000)T(二)原点位移反幂法的MATLAB主程序2用原点位移反幂法计算矩阵A的特征值和对应的特征向量的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('请注意迭代次数k,特征值的近似值lambda,对应的特征向量的近似向量Vk,相邻两次迭代的误差Wc如下：')elsedisp(-请注意迭代次数k已经达到最大迭代次数maxi,特征值的近似值lambda,对应的特征向量的近似向量Vk,相邻两次迭代的误差Wc如下：-)end[V,D]=eig(A,'nobalance'),Vk=U;k=k-1;Wc;lambdan=jlamb+1/mk;例5.3.4用原点位移反幕法的迭代公式(5.27)，计算例题5.3.3,并且将这两个例题的计算结果进行比较.再用两种原点位移反幕法的MATLAB主程序，求「0.99999999999997对应的特征向量.解(1)计算特征值X1^^1=1.001对应特征向量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)运行后屏幕显示结果请注意迭代次数k,特征值的近似值lambda,对应的特征向量的近似向量Vk,相邻两次迭代的误差Wc如下：k=lambda= Wc=5 1.00200000000138 1.376344154436924e-006Vk’= -0.50000000000000 -0.50000000000000-1.00000000000000同理可得，另外与两个特征值对应的特征向量.⑵再用两种原点位移反幂法的MATLAB主程序，求%=0.99999999999997对应的特征向量.输入MATLAB程序>>A=[011-5;-217-7;-426-10];V0=[1,1,1]';[k,lambda,Vk,Wc]=ydwyfmf(A,V0,0.99999999999997,0.001,100)运行后屏幕显示结果请注意：因为A-aE的各阶主子式都不等于零，所以A-aE能进行LU分解.A-aE的秩R(A-aE)和各阶顺序主子式值hl、迭代次数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)运行后屏幕显示结果请注意迭代次数k,特征值的近似值lambda,对应的特征向量的近似向量Vk,相邻两次迭代的误差Wc如下：Vk= k=0.50000000000000 30.50000000000000 lambda=1.00000000000000 1.00000000000000Wc=5.412337245047640e-0165.4雅可比(Jacobi)方法及其MATLAB程序5.4.3雅可比方法的MATLAB程序用雅可比方法计算对称矩阵A的特征值和对应的特征向量的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(-请注意迭代次数k已经达到最大迭代次数max1，迭代次数k,对称矩阵Bk,以特征向量为列向量的矩阵V,特征值为对角元的对角矩阵D如下：')elsedisp(-请注意迭代次数k,对称矩阵Bk,以特征向量为列向量的矩阵V,特征值为对角元的对角矩阵D如下：')endWc;k=k;V=Vk;Bk=B2;D=diag(diag(Bk));[V1,D1]=eig(A,'nobalance')5.常微分方程数值解法用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')得到的结果是解析解没有找到。用数值解。在Matlab下输入：edit，然后将下面两行百分号之间的内容，复制进去，保存%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%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.复化梯形公式用复化梯形公式求解sinx积分，积分区间0—pi建立Trapezoid.m文件function[I,step]=Trapezoid(f,a,b,eps)%f为函数，a为积分上限，b为积分下限，eps为积分精度，step为划分区间个数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次的复化梯形公式积分x=a+h*i;%i=0和n-1时，分别代表积分区间的左右端点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) %程序默认精度为1e-4输出：q=1.9985;s=33若输入[q,s]=Trapezoid(‘sin(x)’，0,6pi,1e输出：q=1.9999;s=157.复化梯形公式求积分(matlab)2007-11-2719:46程序function[I,n,Ichain]=computT(a,b,errorBound,dNum)%复化梯形公式求积分%调用格式：[I,n,Ichain]=computT(a,b,errorBound,dNum)%输入4：%a:积分下限b：积分上限errorBound:输出结果的精度dNum:区间初始分点数%输出3：% I:积分近似值n：最终区间分点数Ichain-迭代过程所有值%被积函数做成函数文件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