《MATLAB源代码.doc》由会员分享，可在线阅读，更多相关《MATLAB源代码.doc（29页珍藏版）》请在金锄头文库上搜索。
1、埃特金插值function f=Atken(x,y,x0)syms t;if(length(x)=length(y) n=length(x);else disp; return;endy1(1:n)=t;for(i=1:n-1) for(j=i+1:n) y1(j)=y(j)*(t-x(i)/(x(j)-x(i)+y(i)*(t-x(j)/(x(i)-x(j); end y=y1 simplify(y1);endif(nargin=3) f=subs(y1(n),t,x0);else simplify(y1(n); f=collect(y1(n); f=vpa(f,6);end特征多项式法fu
2、nction l=Chapoly(A)%特征多项式法求矩阵特征值%已知矩阵：A%求得的矩阵特征值：lsyms t;N=size(A);n=N(1,1);y=det(A-t*eye (n,n);l=solve(y);l=vpa(1,5); %function f = Chebyshev(y,k,x0)syms t;T(1:k+1) = t;T(1) = 1;T(2) = t;c(1:k+1) = 0.0;c(1)=int(subs(y,findsym(sym(y),sym(t)*T(1)/sqrt(1-t2),t,-1,1)/pi;c(2)=2*int(subs(y,findsym(sym(y)
3、,sym(t)*T(2)/sqrt(1-t2),t,-1,1)/pi;f = c(1)+c(2)*t;for i=3:k+1T(i) = 2*t*T(i-1)-T(i-2);c(i) = 2*int(subs(y,findsym(sym(y),sym(t)*T(i)/sqrt(1-t2),t,-1,1)/2;f = f + c(i)*T(i);f = vpa(f,6);if(i=k+1)if(nargin = 3)f = subs(f,t,x0);elsef = vpa(f,6);endendend弦割法求解非线性方程function x k t=ChordsecantToEquation(f
4、,x0,x1,eps)%弦割法求解非线性方程%x k t=ChordsecantToEquation(f,x0,x1,eps)%x:近似解%k:迭代次数%t:运算时间%f:原函数，定义为内联函数%x0,x1:初始值%eps:误差限%应用举例：%f=inline(x3+4*x2-10);%x=ChordsecantToEquation(f,1,2,0.5e-6)%x k=ChordsecantToEquation(f,1,2,0.5e-6)%x k t=ChordsecantToEquation(f,1,2,0.5e-6)%函数的最后一个参数也可以不写，默认情况下，eps=0.5e-6%x k
5、t=ChordsecantToEquation(f,1,2)if nargin=3eps=0.5e-6;endtic;k=0;while 1x=x1-f(x1)*(x1-x0)./(f(x1)-f(x0);k=k+1;if abs(x-x1) 30break;endx0=x1;x1=x;endt=toc;if k = 30disp(迭代次数太多。);x=0;t=0;end复合求积公式function I,step = CombineTraprl(f,a,b,eps)%f 被积函数%a,b 积分上下限%eps 精度%I 积分结果%step 积分的子区间数if(nargin =3)eps=1.0e
6、-4;endn=1;h=(b-a)/2;I1=0;I2=(subs(sym(f),findsym(sym(f),a)+subs(sym(f),findsym(sym(f),b)/h;while abs(I2-I1)epsn=n+1; h=(b-a)/n;I1=I2;I2=0;for i=0:n-1x=a+h*i;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;function f = DCS(x,y,x0)syms t;if(length(x)
7、= length(y)n = length(x);c(1:n) = 0.0;elsedisp(x和y的维数不相等！);return;endc(1) = y(1);for(i=1:n-1) for(j=i+1:n)y1(j) = (x(j)-x(i)/(y(j)-y(i);endc(i+1) = y1(i+1); y = y1; endf = c(n);for(i=1:n-1)f = c(n-i) + (t-x(n-i)/f;f = vpa(f,6);if(i=n-1)if(nargin = 3)f = subs(f,t,x0);elsef = vpa(f,6);endendend;用欧拉法求一
8、阶常微分方程的数值解 function y = DEEuler(f, h,a,b,y0,varvec)format long;N = (b-a)/h;y = zeros(N+1,1);x = a:h:b;y(1) = y0;for i=2:N+1y(i) = y(i-1)+h*Funval(f,varvec,x(i-1), y(i-1);endformat short; 用隐式欧拉法求一阶常微分方程的数值解 function y = DEimpEuler(f, h,a,b,y0,varvec)format long;N = (b-a)/h;y = zeros(N+1,1);y(1) = y0;
9、x = a:h:b;var = findsym(f);for i=2:N+1fx = Funval(f,var(1),x(i);gx = y(i-1)+h*fx - varvec(2);y(i) = NewtonRoot(gx,-10,10,eps);endformat short;龙格库塔法function y = DELGKT2_suen(f, h,a,b,y0,varvec)format long;N = uint16(b-a)/h);y = zeros(N+1,1);y(1) = y0;x = a:h:b;var = findsym(f);for i=2:N+1 K1 = Funval
10、(f,varvec,x(i-1) y(i-1); K2 = Funval(f,varvec,x(i-1)+2*h/3 y(i-1)+K1*2*h/3); y(i) = y(i-1)+h*(K1+3*K2)/4;endformat short; .用改进欧拉法求一阶常微分方程的数值解 function y = DEModifEuler(f, h,a,b,y0,varvec)format long;N = (b-a)/h;y = zeros(N+1,1);y(1) = y0;x = a:h:b;var = findsym(f);for i=2:N+1v1 = Funval(f,varvec,x(i
11、-1) y(i-1);t = y(i-1) + h*v1;v2 = Funval(f,varvec,x(i) t);y(i) = y(i-1)+h*(v1+v2)/2;endformat short; 一阶微分方程Function dy=myfun_1(x,y) %定义输入,输出变量和函数文件名 dy=zeros(2,1); %明确dy的维数,用微分方程组时不可缺省 dy(1)=y(2); %dy(m)表示y的m阶导数;y(n)表示y的第n列 dy(2)=2*y(2)-2*y(1)+5*exp(2*x)*sin(x); %与方程组中第二个微分方程对应 基本Guass消去法求解线性方程组func
12、tion x=EqtsBasicGuass(A,b)%基本Guass消去法求解线性方程组Ax=b%x=EqtsBasicGuass(A,b)%x:解向量，列向量%A:线性方程组的矩阵%b:列向量%检查输入参数if size(A,1) = size(b,1)disp(输入参数有误！);x= ;return;end%(A|b)A=A b;%消去过程n=size(A,1);l=zeros(n);for k=1:n-1for i=k+1:nl(i,k)=A(i,k)/A(k,k);endfor i=k+1:nfor j=k+1:n+1A(i,j)=A(i,j)-l(i,k)*A(k,j);endfor
13、 j=1:kA(i,j)=0;endendend%回代过程x=zeros(n,1);x(n)=A(n,n+1)/A(n,n);for i=n-1:-1:1y=0;for j=i+1:ny=y+A(i,j)*x(j);endx(i)=(A(i,n+1)-y)/A(i,i);endreturn;追赶法求解三对角线性方程组function x=EqtsForwardAndBackward(L,D,U,b)%追赶法求解三对角线性方程组Ax=b%x=EqtsForwardAndBackward(L,D,U,b)%x:三对角线性方程组的解%L:三对角矩阵的下对角线，行向量%D:三对角矩阵的对角线，行向量%U:三对角矩阵的上对角线，行向量%b:线性方程组Ax=b中的b，列向量%应用举例:%L=-1 -2 -3;D=2 3 4 5;U=-1 -2 -3;b=6 1 -2 1;%x=EqtsForwardAndBackward(L,D,U,b)%检查参数的输入是否正确n=length(D);m=length(b);n1=length(L);n2=length(U);if n-n1 = 1 | n-n2 = 1 | n = mdisp(输入参数有误！)x= ;return;end%追的过程for i=2:nL(i-1)
