关键函数.docx
- 文档编号:10565918
- 上传时间:2023-05-26
- 格式:DOCX
- 页数:22
- 大小:20.65KB
关键函数.docx
《关键函数.docx》由会员分享,可在线阅读,更多相关《关键函数.docx(22页珍藏版)》请在冰点文库上搜索。
关键函数
1,逐次超松弛迭代法:
functiony=SORit(w,A,b,x,it_limit)
fprintf('逐次超松弛迭代法\n');
[n,n]=size(A);
D=diag(diag(A));
L=-tril(A,-1);
U=-triu(A,1);
%迭代矩阵
B=inv(D-w*L)*((1-w)*D+w*U)
b=w*inv(D-w*L)*b
it=0;
while1
it=it+1;
if(it>it_limit)
fprintf('迭代次数超界\n');
break
end
y=B*x+b;
if(abs(y-x)<1e-3)
fprintf('迭代次数=%i\n',it);
break;
else
x=y;
end
end
it
v=[1,1,1]
v=
111
>>diag(v)
ans=
100
010
001
>>diag(diag(v))
ans
1
1
1
用法:
tril(X),其中X表示一个矩阵。
功能:
tril函数是matlab中提取矩阵下三角矩阵的函数。
tril(X)产生X矩阵的下三角矩阵,其余元素补0。
相关函数:
triu函数提取矩阵的上三角矩阵。
第二章例子:
二分法:
functionbisec2(f_name,a,c)
%f_name:
需求根方程函数名
%a:
左端点
%c:
右端点
%%tolerance:
精度要求
%it_limit:
最大迭代次数控制
fprintf('\n二分法求解方程%s=0\n',f_name);
tolerance=1e-6;
it_limit=30;
it=0;
Y_a=feval(f_name,a);
Y_c=feval(f_name,c);
if(Y_a*Y_c>0)
fprintf('不满足条件f(a)f(c)<=0\n');
else
while1
it=it+1;
b=(a+c)/2;
Y_b=feval(f_name,b);
if(abs(c-a)<=tolerance)
fprintf('根=%14.8f\n迭代次数=%i\n',b,it);
break%达到精度要求
end
if(it>it_limit)
fprintf('迭代次数超界\n');break
end
if(Y_a*Y_b<=0)
c=b;Y_c=Y_b;
else
a=b;Y_a=Y_b;
end
end
end
简单迭代法:
functiondiedai(x0)
%x0表示初值
%用迭代方程解
%%tolerance:
精度要求
%it_limit:
最大迭代次数控制
fprintf('\n用迭代方式求解xk+1=(2-exp(xk))/10');
tolerance=0.5*10.^-3;
it_limit=30;
x1=x0;
it=0;
while1
it=it+1;
x2=(2-exp(x1))/10;
fprintf('\nRoot=%14.8f\n迭代次数=%i\n',x2,it);
if(abs(x2-x1) break%达到精度要求 end if(it>it_limit) fprintf('迭代次数超界\n');break end x1=x2; end functiony=simpleit(A,b,x,it_limit) fprintf('简单迭代法\n'); [n,n]=size(A); B=eye(n)-A it=0; while1 it=it+1; if(it>it_limit) fprintf('迭代次数超界\n'); break end y=B*x+b; if(norm(y-x)<1e-5) break; else x=y; end end it 牛顿迭代法: functionnewt2(f_name,f1_name,x0) %f_name: 所求函数 %f1_name: 所求函数的导数 %x0: 初值 %%tolerance: 精度要求 %it_limit: 最大迭代次数控制 fprintf('\n用迭代方式求解xk+1=xk-f(x)/f;(x)'); tolerance=0.5*10.^-3; it_limit=30; x1=x0; it=1; while1 y=feval(f_name,x1); y1=feval(f1_name,x1); x2=x1-y/y1; fprintf('\n迭代次数: it=%3.0fRoot=%12.6f\n',it,x2); if(abs(x2-x1) break end if(it>it_limit) fprintf('迭代次数超界\n');break end x1=x2; it=it+1; end 第三章例子: 高斯消元法: functionb=Gauss2(A,b) %求解AX=b,结果保存在b中(这是老师的思想,觉得这样是一个好方法就直接用了) %A为系数矩阵,b为最右端的列向量 %i=k+1: n %A(k,k)=a(kk) [n,n]=size(A); %[m,n]=size(X)returnsthesizeofmatrixXinseparatevariablesmandn %: size函数直接求出系数矩阵A的行列大小 fprintf('高斯消元法求Ax=b的解\n'); %%第一步消元运算 fork=1: n-1 %一共要进行n-1次消元就可以转换为上三角矩阵 A(k+1: n,k)=A(k+1: n,k)/A(k,k); %Mik=A(k+1: n,k)求出消去未知元的系数 %注意i=k+1: n在k=某一个值的时候将做一次循环,求出i为不同值的时候Mik的值,下同 A(k+1: n,k+1: n)=A(k+1: n,k+1: n)-A(k+1: n,k)*A(k,k+1: n) %消元并刷新值 b(k+1: n)=b(k+1: n)-A(k+1: n,k)*b(k); %刷新的是b中的值,因为要得到相同的解,所以右端列矩阵要做相同的变换 end %%第二步回代求解 b(n)=b(n)/A(n,n); %求出最后一个解 fork=n-1: -1: 1 b(k)=(b(k)-A(k,k+1: n)*b(k+1: n))/A(k,k); %从后往前依次求出解 End 列主元消元法 functionb=liezhuxiaoyuan(A,b) %该函数是用列主消元法解Ax=b的解, %其中A是系数矩阵,b是右端列向量,而解也最后保存在b中 fprintf('列主消元法\n'); %求系数矩阵的阶保存在n中 [n,n]=size(A); piv=1: n; %%: 经过n-1步就完成变换 %%选主元: 选择ik使得|a(ik)|=max fork=1: n-1 %行保存在maxv中,列保存在r中 [maxv,r]=max(abs(A(k: n,k))); ifmaxv<1e-6 %违背了列主元消元法的条件与初衷 fprintf('主元过小,停止计算\n'); break; return end %%换行消元操作 q=r+k-1; piv([kq])=piv([qk]); A([kq],: )=A([qk],: ); b([kq])=b([qk]); A(k+1: n,k)=A(k+1: n,k)/A(k,k); A(k+1: n,k+1: n)=A(k+1: n,k+1: n)-A(k+1: n,k)*A(k,k+1: n); b(k+1: n)=b(k+1: n)-A(k+1: n,k)*b(k); A end %%逆序求解存放在b中 b(n)=b(n)/A(n,n); fork=n-1: -1: 1 b(k)=(b(k)-A(k,k+1: n)*b(k+1: n))/A(k,k); end 雅克比迭代法 functiony=MyJacobiit(A,b,x,it_limit) %A为系数矩阵。 %b为常数项组成的向量。 %x为初始向量X(0)。 %it_limit为用户设置的迭代次数的上限值。 fprintf('雅可比迭代法\n'); [n,n]=size(A); D=diag(diag(A)); B=inv(D)*(D-A)%B为化简后的便于迭代的等价方程组的系数矩阵。 b=inv(D)*b%b为化简后的便于迭代的等价方程组的常数项向量。 it=0; while1 it=it+1; if(it>it_limit) fprintf('迭代次数超界\n'); break end y=B*x+b; if(norm(y-x)<1e-4)%近似于条件||X(k+1)-X(k)||的2范数<=0.0001。 break; else k=it%k为当前迭代次数。 x=y end end all=it%all为实际总迭代次数。 functiony=Jacobiit(A,b,x,it_limit) fprintf('雅可比迭代法\n'); %x为初始向量 %求矩阵的阶数 [n,n]=size(A); %diag表示取对角元 %diag(diag(A))表示取出的对角元为对角的对角矩阵(除对角外都为0) D=diag(diag(A)); %inv用于求逆矩阵 %A=D-L-U %L+U=D-A %雅克比迭代矩阵Bj=inv(D)(L+U)fj=inv(D)*b B=inv(D)*(D-A) %雅克比矩阵 b=inv(D)*b %fj; it=0; while1 it=it+1; if(it>it_limit) fprintf('迭代次数超界\n'); break end y=B*x+b; %norm(y-x,inf) if(norm(y-x,inf)<1e-3) fprintf('迭代次数=%i\n',it); break; else x=y; end end it 高斯-赛德尔迭代法 functiony=MyGaussSeidelit(A,b,x,it_limit) %A为系数矩阵。 %b为常数项组成的向量。 %x为初始向量X(0)。 %it_limit为用户设置的迭代次数的上限值。 fprintf('高斯-赛德尔迭代法\n'); [n,n]=size(A); DL=tril(A); B=inv(DL)*(DL-A)%B为化简后的等价方程组的系数矩阵(并非便于迭代的)。 b=inv(DL)*b%b为化简后的等价方程组的常数项向量(并非便于迭代的)。 %%将B与b带入方程组,化简约项后才可得高斯-赛德尔迭代法的对应的便于迭代的方程组。 it=0; while1 it=it+1; if(it>it_limit) fprintf('迭代次数超界\n'); break end y=B*x+b; if(norm(y-x)<1e-4)%近似于条件||X(k+1)-X(k)||的2范数<=0.0001。 break; else k=it%k为当前迭代次数。 x=y end end all=it%all为实际总迭代次数。 逐次超松弛 functiony=MySORit(w,A,b,x,it_limit) %w为松弛因子。 %A为系数矩阵。 %b为常数项组成的向量。 %x为初始向量X(0)。 %it_limit为用户设置的迭代次数的上限值。 fprintf('逐次超松弛迭代法\n'); [n,n]=size(A); D=diag(diag(A)); L=-tril(A,-1); U=-triu(A,1); B=inv(D-w*L)*((1-w)*D+w*U) b=w*inv(D-w*L)*b it=0; while1 it=it+1; if(it>it_limit) fprintf('迭代次数超界\n'); break end y=B*x+b; if(norm(y-x)<1e-4)%近似于条件||X(k+1)-X(k)||的2范数<=0.0001。 break; else %k=it%k为当前迭代次数。 x=y;%此时不输出每次迭代后向量X的值,因为当w接近2时,迭代次数可能会很多。 end end all=it%all为实际总迭代次数。 第四章例子: 三对角追赶法: %求解三对角线方程组的追赶法(n阶) %韩臻,2003年5月16日编制 %参见教材第51-53页。 %输入参数: 系数矩阵A的阶数n,三对角线上的元素数组a,b,c,和右端项f, %它们都是一维向量数据(长为n),其中b是主对角线,a是下对角线; %返回参数: f是求得的解。 functionf=SDJ1(n,a,b,c,f)%定义追赶法函数 eps=0.1e-15;%除数为0的判别阈值 %以下计算alpha和beta ifabs(b (1)) disp('除数为0,停止计算'); return else c (1)=c (1)/b (1); end fori=2: n-1 b(i)=b(i)-a(i)*c(i-1); ifabs(b(i)) disp('除数为0,停止计算'); return else c(i)=c(i)/b(i); end end b(n)=b(n)-a(n)*c(n-1); %以下求解Ly=f,y的值存入f f (1)=f (1)/b (1); fori=2: n f(i)=(f(i)-a(i)*f(i-1))/b(i); end %以下求解Ux=y,x的值存入f fori=n-1: -1: 1 f(i)=f(i)-c(i)*f(i+1); end return %结束 。 拉格朗日插值法: %Lagrange插值多项式(n次) %韩臻,2003年4月14日编制 %参见《科学计算引论(第二版)》(Nakamura著,梁恒等译)第145页) %输入参数: 插值数据(x,f),其中x和f是一维向量数据; %待计算的点的横坐标数组xi; %返回参数: fi是用Lagrange插值多项式P(x)求得的P(xi)值. functionfi=Lagran_(x,f,xi) fi=zeros(size(xi)); np1=length(f); fori=1: np1 z=ones(size(xi)); forj=1: np1 ifi~=j,z=z.*(xi-x(j))/(x(i)-x(j));end end fi=fi+z*f(i); end return %结束。 %%高次(n)多项式拟合曲线例子、、最小二乘法 x=[0.1,0.4,0.5,0.7,0.7,0.9]; y=[0.61,0.92,0.99,1.52,1.47,2.03]; %% n=3; c=polyfit(x,y,n) %% x=x'; y=y'; L=length(x); %A(: n+1)=ones(L,1); %fork=n: -1: 1 %A(: k)=x.*A(: k+1); %end %cc=A\y %c3=(A'*A)\(A'*y) xx=0: 0.01: 1; yy=polyval(c,xx); plot(x,y,'r*',xx,yy,'-m') %% y1=polyval(c,x); fprintf('平方误差=%f\n',norm(y-y1)); fprintf('最大偏差=%f',norm(y-y1,inf)); %%直线拟合例子 %p.95例5 x=[1,2,4,6,8,10]; y=[1.8,3.7,8.2,12.0,15.8,20.2]; plot(x,y,'ro') axis([012,-525]) pause holdon c=polyfit(x,y,1) xx=0: 0.1: 11; yy=c (1)*xx+c (2); plot(xx,yy) %% x=x'; y=y'; A(: 1)=x A(: 2)=1 cc=A\y c3=(A'*A)\(A'*y) AA=A'*A y1=A'*y %%一般函数线性组合拟合例子最小二乘法 %y=c1+c2*x+c3*sin(x)+c4*exp(x) x=[0.1,0.4,0.5,0.7,0.7,0.9]; y=[0.61,0.92,0.99,1.52,1.47,2.03]; %% x=x'; y=y'; L=length(x); A(: 1)=ones(L,1); A(: 2)=x; A(: 3)=sin(x); A(: 4)=exp(x); cc=A\y c3=(A'*A)\(A'*y) xx=0: 0.01: 1; yy=cc (1)+cc (2)*xx+cc(3)*sin(xx)+cc(4)*exp(xx); holdon plot(x,y,'r*',xx,yy) %axis([01,03]) %% y1=cc (1)+cc (2)*x+cc(3)*sin(x)+cc(4)*exp(x); fprintf('平方误差=%f\n',norm(y-y1)); fprintf('最大偏差=%f',norm(y-y1,inf)); 第五章 %%复合梯形数值积分公式-满足精度eps要求 %%2011年修订 %%n是初始取点数,一般可以取大一些,以防止假收敛。 functionTn=Ttest(fun,a,b,eps,n) h=(b-a)/n; x=a: h: b; y=feval(fun,x); T0=h*(sum(y)-0.5*(y (1)+y(n+1))); while1 n=2*n; h=0.5*h; x2=a+h: 2*h: b-h; y2=feval(fun,x2); Tn=0.5*T0+h*sum(y2); ifabs(Tn-T0) break; else T0=Tn; end end n Tn %% functionsn=simpson(fun,a,b,n) %复合辛普森公式 %a为积分下限,b为积分上限,n为将[a,b]分为n等分 h=(b-a)/n; %% x1=a: h: b %[y1,y2,...]=feval(fhandle,x1,...,xn) %evaluatesthefunctionhandle,fhandle,usingargumentsx1throughxn. y=feval(fun,x1); %% x2=a+0.5*h: h: b-0.5*h; y2=feval(fun,x2); %%求sn sn=h*(y (1)+2*sum(y(2: n))+4*sum(y2)+y(n+1))/6; %% %%复合辛普森数值积分公式 %%2010年修订 functionSn=simpson_1(fun,a,b,n) h=(b-a)/n; x=a: h: b; y=feval(fun,x); %复合辛普森公式Sn x2=a+0.5*h: h: b-0.5*h; y2=feval(fun,x2); Sn=h*(y (1)+4*sum(y2)+2*sum(y(2: n))+y(n+1))/6; %% %%复合辛普森数值积分公式,满足精度eps %%2011年修订 functionSn=simpson_2(fun,a,b,eps) n=8; S0=simpson_1(fun,a,b,n); while1 n=2*n; S1=simpson_1(fun,a,b,n); ifabs(S1-S0) break; end S0=S1; end n S1 %% %%复合低阶牛顿-科茨数值积分公式试验1 %%n相同的复合梯形、辛普森、科茨公式 %%2010年修订 functionTSCtest1(fun,a,b,n) h=(b-a)/n; x=a: h: b; y=feval(fun,x); %% %复合梯形公式Tn Tn=0.5*h*(y (1)+2*sum(y(2: n))+y(n+1)); %% %复合辛普森公式Sn x2=a+0.5*h: h: b-0.5*h; y2=feval(fun,x2); Sn=h*(y (1)+4*sum(y2)+2*sum(y(2: n))+y(n+1))/6; %% %复合科茨公式Cn x4=a+0.25*h: h: b-0.75*h; x3=a+0.75*h: h: b-0.25*h; y4=feval(fun,x4); y3=feval(fun,x3); Cn=h*(7*y (1)+32*sum(y4)+12*sum(y2)+32*sum(y3)+14*sum(y(2: n))+7*y(n+1))/90; %% fprintf('T%d=%20.16f用了%d个结点\n',n,Tn,n+1); fprintf('S%d=%20.16f用了%d个结点\n',n,Sn,2*n+1); fprintf('C%d=%20.16f用了%d个结点\n',n,Cn,4*n+1);
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 关键 函数
![提示](https://static.bingdoc.com/images/bang_tan.gif)