matlab程序运行修改版

1、拉格朗日插值公式：function y=lage(X,Y,x)% 拉格朗日插值函数n=length(X);y=zeros(size(x);for k=1:n w=ones(size(x); for j=1:k-1 k+1:n w=(x-X(j)./(X(k)-X(j).*w; end y=y+w*Y(k);end调用方法：在matlab中新建一个script窗口，把上述程序保存成（lage.m）文件,然后在命令窗口调用y=lage(A,B,x),其中A，B均为列向量，分别代表插值点的横坐标和纵坐标，x为待求点横坐标例如:现有x0=1,x1=2,x2=4,y0=3,y1=2,y2=7;要求x=1.5处的估计值。过程：则A=1 2 3,B=3 2 7，x=1.5在命令窗口输入如下命令：A=1 2 3;B=3 2 7;x=1.5;y=lage(A,B,x)(注意：这个直接占到命令窗口会出现符号问题，请自行修改)y即为所求结果2、牛顿插值公式：function f=niudun(X,Y,x)n=length(X);A=zeros(n,n);%A用于存储均差表（空余位置设为0）A(:,1)=Y;t=A(1,1); for j=2:n for i=j:n l=1; A(i,j)=(A(i,j-1)-A(i-1,j-1)./(X(i)-X(i-j+1);%算出均差表 if i=j%只取每列非0的第一个数 for k=1:j-1 l=l.*(x-X(k);%通项计算 end t=t+l.*A(i,j);%各项累加求结果 end endendf=t;在matlab中新建一个script窗口，把上述程序保存成(niudun.m)调用方法：y=niudun(A,B,x) A,B,x的定义与拉格朗日函数相同3、 雅克比迭代法：function y,n=yacobi(a,d)format long;x=0;0;0; %stop=1.0e-4 ; %L=-tril(a,-1);U=-triu(a,1); A=(diag(diag(a);D=inv(A);X=D*(L+U)*x+D*d; % Jn=1;while norm(X-x,inf)=stop % x=X; X=D*(U+L)*x+D*d;n=n+1;endy=X;或：function yacobia=3 0 -2;0 8 1;-1 1 6 ;d=3;2;2;x=0;0;0; %初始向量stop=1.0e-4 ; %迭代的精度L=-tril(a,-1);U=-triu(a,1); A=(diag(diag(a);D=inv(A);X=D*(L+U)*x+D*d; % J迭代公式n=1;while norm(X-x,inf)=stop % 时迭代中止否则继续x=X; X=D*(U+L)*x+D*d;n=n+1;endXn在matlab中新建一个script窗口，在这个窗口中（即.m文件中）直接运行（run）上述程序，若题目要求Ax=b,则将a矩阵内容换为A，d内容换为b。例如：数值分析书上209页的第一题过程：将上述程序二、三行改为a=5 2 1;-1 4 2;2 -3 10; d=-12;20;3其余不变，即可得到迭代次数与结果。4、 高斯赛德尔迭代法：function gaosisaia=3 0 -2;0 8 1;-1 1 6 ;d=3;2;2;x=0;0;0; %stop=1.0e-4; %L=-tril(a,-1);U=-triu(a,1); A=(diag(diag(a);D=inv(A-L);X=D*U*x+D*d; % Jn=1;while norm(X-x,inf)=stop % x=X; X=D*U*x+D*d;n=n+1;endXn调用方法：在matlab中新建一个script窗口，在这个窗口中（即.m文件中）直接运行（run）上述程序，若题目要求Ax=b,则将a矩阵内容换为A，d内容换为b。（同3）5、高斯消去法：function gaosia=0 1 2 3;9 11 23 34;62.5 23.4 15.5 17.2;192.01 124 25.1 59.3;b=1;1;1;1;n=length(b);m=zeros(n,n);for k=1:n-1 m=k+1; while(a(k,k)=0&m=n) c=a(k,:);a(k,:)=a(k+1,:);a(k+1,:)=c;d=b(k,:);b(k,:)=b(k+1,:);b(k+1,:)=d;endfor i=k+1:n m(i,k)=a(i,k)/a(k,k); b(i)=b(i)-m(i,k)*b(k); for j=k:n a(i,j)=a(i,j)-m(i,k)*a(k,j); endendendx(n)=b(n)/a(n,n);for k=n-1:-1:1 s=0; for j=k+1:n s=s+a(k,j)*x(j); end x(k)=(b(k)-s)/a(k,k);endx=x调用方法：在matlab中新建一个script窗口，在这个窗口中（即.m文件中）直接运行（run）上述程序，若题目要求Ax=b,则将a矩阵内容换为A，d内容换为b。（同3、4）6、高斯列主元消去法：function liezhua=0 1 2 3;9 11 23 34;62.5 23.4 15.5 17.2;192.01 124 25.1 59.3 ;b=1;1;1;1;n=length(b);det=1;for k=1:n-1 c,u=max(abs(a(k:n,k); %取出最大值和位置 ik=u+k-1; %最大值在举证中的位置 if c=0; %如果最大值为0，直接结束循环 det=0; break; end if ik=k bk=b(k);b(k)=b(ik);b(ik)=bk; ck=a(k,:);a(k,:)=a(ik,:);a(ik,:)=ck; end for i=k+1:n %这里开始和高斯消去一样 m(i,k)=a(i,k)/a(k,k); b(i)=b(i)-m(i,k)*b(k); for j=k:n a(i,j)=a(i,j)-m(i,k)*a(k,j); end endendif det=0 for k=1:n det=det*a(k,k); endendx(n)=b(n)/a(n,n);for k=n-1:-1:1 s=0; for j=k+1:n s=s+a(k,j)*x(j); end x(k)=(b(k)-s)/a(k,k);endx=x调用方法：在matlab中新建一个script窗口，在这个窗口中（即.m文件中）直接运行（run）上述程序，若题目要求Ax=b,则将a矩阵内容换为A，d内容换为b。（同3、4、5）7、牛顿迭代法：function x,n=newton(x0,e,N)%x=newton(e,et,n)n=0; while fdiff(x0)=0|n=N x1=x0-f(x0)/fdiff(x0); if abs(x1-x0)e x=x1; return end x0=x1; n=n+1;end endfunction y=f(x) y=x2-3*x+2-exp(x); endfunction y=fdiff(x0)x=sym(x);%sym是将字符或者数字转换为字符y=diff(x2-3*x+2-exp(x);y=subs(y,x0);%subs是赋值函数,用数值替代符号变量替换函数end说明：在matlab中新建一个script窗口，写入上述程序，上述文件中含有三个fuction，其中后两个中的函数（x2-3*x+2-exp(x)）即是要求求解的f(x)=0中的f(x)，实际操作时应先把这两个函数改为题目要求的函数，然后将其保存为newton.m文件。调用方法：在命令窗口输入x n=newton(0,10-4,1000)(其中0可改为任意实数，但最好不要离精确解太远；10-4为迭代精度，依据需要进行修改；1000为最大迭代次数，一般情况下可不变)8、二分法：（注意：a,b必须是f(x)上异号的两点，否则极有可能运行错误）function x,n=erfen(a,b,e,N)n=0;x1=a;x2=b;if f(a)=0x0=a;else if f(b)=0 x0=b; else while n=N x0=(x1+x2)/2; if abs(f(x0)e x=x0; return else if f(x0)*f(x1)0.000005) n=n*2; a=comsimpson(fun,0,1,n); b=comsimpson(fun,0,1,2*n);endnaendfunction y=comsimpson(fun,a,b,n) z1=feval (fun,a)+ feval (fun,b);m=n/2;h=(b-a)/(2*m);x=a;z2=0; z3=0; x2=0; x3=0;for k=2:2:2*m x2=x+k*h; z2= z2+2*feval (fun,x2); endfor k=1:2:2*m x3=x+k*h; z3= z3+4*feval (fun,x3); endy=(z1+z2+z3)*h/3;endfunction y=fun(x)y=1/(x2+3);end将上述程序保存为一个a.m文件，把最后一个函数换掉即可 梯形公式：function an=4;a=tixing(fun,0,1,n);b=tixing(fun,0,1,2*n);while(abs(b-a)0.000005) n=n*2; a= tixing (fun,0,1,n); b= tixing (fun,0,1,2*n);endnaendfunction y=tixing(fun,a,b,n) z1=feval (fun,a)+ feval (fun,b);h=(b-a)/n;x=a;z2=0; x2=0; for k=1:n-1 x2=x+k*h; z2= z2+2*feval (fun,x2); endy=(z1+z2)*h/2;endfunction y=fun(x)y=1/(x2+3);end方法二：完成这两个公式共三个步骤：（1）构造待积分函数：function y=fun(x)y=1/(x2+3);调用方法：在matlab中新建一个script窗口，把上述程序保存成fun.m文件，关闭。（2）构造复合梯形公式或复合辛普森公式：复合梯形公式如下：function y=tixing(fun,a,b,n) z1=feval (fun,a)+ feval (fun,b);h=(b-a)/n;x=a;z2=0; x2=0; for k=1:n-1 x2=x+k*h; z2= z2+2*feval (fun,x2); endy=(z1+z2)*h/2;调用方法：在matlab中新建一个script窗口，把上述程序保存成保存为tixing.m文件，关闭。复合辛普森公式如下：function y=comsimpson(fun,a,b,n) z1=feval (fun,a)+ feval (fun,b);m=n/2;h=(b-a)/(2*m);x=a;z2=0; z3=0; x2=0; x3=0;for k=2:2:2*m x2=x+k*h; z2= z2+2*feval (fun,x2); endfor k=1:2:2*m x3=x+k*h; z3= z3+4*feval (fun,x3); endy=(z1+z2+z3)*h/3;在matlab中新建一个script窗口，把上述程序保存成comsimpson.m文件，关闭。（3） 迭代求出近似解：n=4;a=comsimpson(fun,0,1,n);b=comsimpson(fun,0,1,2*n);while(abs(b-a)0.000005) n=n*2; a=comsimpson(fun,0,1,n); b=comsimpson(fun,0,1,2*n);endna（梯形公式即将comsimpson改为tixing,0.000005为精度，根据题目修改）调用方法：在matlab中新建一个script窗口，写入上述程序，并运行，在命令窗口可得到结果 9