数值分析编程总结.docx
- 文档编号:4387056
- 上传时间:2023-05-07
- 格式:DOCX
- 页数:22
- 大小:26.26KB
数值分析编程总结.docx
《数值分析编程总结.docx》由会员分享,可在线阅读,更多相关《数值分析编程总结.docx(22页珍藏版)》请在冰点文库上搜索。
数值分析编程总结
1.LU分解:
[LU]=lu(A);
2.追赶法
functionx=zhuiganfa(A,b)
[n,n]=size(A);
fori=1:
n
if(i==1)
l(i)=A(i,i);
y(i)=b(i)/l(i);
else
l(i)=A(i,i)-A(i,i-1)*u(i-1);
y(i)=(b(i)-y(i-1)*A(i,i-1))/l(i);
end
if(i u(i)=A(i,i+1)/l(i); end end x(n)=y(n) forj=n-1: -1: 1 x(j)=y(j)-u(j)*x(j+1); end 数值试验: n=101; a=1210………………………………..….0 1121………………………………..…0 01121………………………..…0 001121……………………0 1 112 关键是如何定义上述矩阵: >>n=101; c1=ones(1,n-1); a1=diag(c1,-1);这个-1说明行位置-1 c2=12*ones(1,n); a2=diag(c2); c3=ones(1,n-1); a3=diag(c3,1); a=a1+a2+a3; 3.拉格朗日插值 functionyh=lage(x,y,xh) n=length(x); m=length(xh); yh=zeros(1,m); c1=ones(n-1,1); c2=ones(1,m); fori=1: n xp=x([1: i-1i+1: n]); yh=yh+y(i)*prod((c1*xh-xp'*c2)./(x(i)-xp'*c2)); end end >>x=[11,12]; >>y=[2,4]; >>xh=[11.75]; >>lage(x,y,xh) ans= 3.5000 4最小二乘法 1.最小二乘的xi和yi为: xi 1953 1964 1982 1990 2000 yi 5.82 6.95 10.08 11.34 12.66 要拟合的函数为: y=a+bx-cxy注意不是多项式 2.编程函数为: functionz=erchen(x,y) x1=ones(5,1); A=[x1,x,-x.*y];注意点乘 z=A\y;注意左除 a=z (1); b=z (2); c=z(3); end 输入: ≻≻x=[19531964198219902000]'; ≻≻y=[5.826.9510.0811.3412.66]'; ≻≻erchen(x,y) ans= 2.9456=a -0.0014=b -0.0005=c 1.最小二乘的xi和yi为: xi 0 0.25 0.5 0.75 1 yi 1 1.284 1.6487 2.1170 2.7183 要拟合的函数为: y=a+bx+cx2是多项式 2.编程函数为: functionz=erchen2(x,y) x1=ones(5,1); A=[x1,x,x.^2]; z=A\y; a=z (1); b=z (2); c=z(3); end 输入: ≻≻x=[00.250.50.751.00]'; ≻≻y=[1.001.2841.64872.11702.7183]'; ≻≻erchen2(x,y) ans= 1.0051 0.8642 0.8437 最小二乘多项式拟合的简单函数方法: ≻≻x=[00.250.50.751.00]'; ≻≻y=[1.001.2841.64872.11702.7183]'; ≻≻P=polyfit(x,y,2)要拟合成4次,则2改成4就可以了 P= 0.84370.86421.0051注意此内置函数输出的结果c,b,a是反的 5复合辛普森公式求解积分 先定义函数: functionv=f(x) v=sin(x);“若定义有除数要点除,分母有0时要特殊定义” end 定义程序: functionI=fsps(f,a,b,n) h=(b-a)/n; x=linspace(a,b,2*n+1); y=feval(f,x); I=(h/6)*(y (1)+2*sum(y(3: 2: 2*n-1))+4*sum(y(2: 2: 2*n))+y(2*n+1)); end >>fsps('f',0,1,4) ans= 0.4597 6.不动点迭代思路 不动点迭代常常有好几个迭代的不动点函数,所以要分别定义这些函数是很困难的,如是乎使用SWITCH内置函数进行切换,叫切换函数. 1.先定义函数后进行编程的方法 先需要定义不动点函数: functionv=f(x) v=x^3-x-1; end 再定义编程: function[it,x]=fixpnt1(f,a,maxit,tol) it=0; x=feval(f,a); whileit<=maxit&abs(x-a)>tol, it=it+1; a=x; x=feval(f,a); end 此函数的调用: >>fixpnt1('f',2,100,1e-5) ans= 1 3.利用切换函数SWITCH的方法(多个不动点迭代函数) function[x,it]=fixpnt(np,a,maxit,tol) switchnp, case1, phi=inline('(3*x+10)^(1/5)'); case2, phi=inline('sin(10*x)+2*cos(x)-3'); case3, phi=inline('3-atan(x)'); case4, phi=inline('-2-1/log(x^2+x+1)'); end it=0; x=phi(a); whileit<=maxit&abs(x-a)>tol, it=it+1; a=x; x=phi(a); end 使用与输入: >>fixpnt(2,1,100,1e-5) ans= -4.2696 7.雅可比迭代 function[xit]=jacobi(A,b,tol) D=diag(diag(A)); L=D-tril(A); U=D-triu(A); x=zeros(size(b)); forit=1: 500 x=D\(b+L*x+U*x); error=norm(b-A*x)/norm(b); if(error break; end end 8.高斯迭代 function[xit]=gaosi(A,b,tol) D=diag(diag(A)); L=D-tril(A); U=D-triu(A); x=zeros(size(b)); forit=1: 500 x=(D-L)\(b+U*x); error=norm(b-A*x)/norm(b); if(error break; end end 9.SOR迭代 function[xit]=SOR(A,b,w,tol) D=diag(diag(A)); L=D-tril(A); U=D-triu(A); x=zeros(size(b)); forit=1: 500 x=(D-w*L)\(w*b+(1-w)*D*x+w*U*x); error=norm(b-A*x)/norm(b); if(error break; end end 10.二分法: 1.先要定义所求的函数: functionv=f(x) v=x^3-x-1; end 2.二分法程序如下: function[x,it]=erfengfa(a,b,f,tol) fa=feval(f,a); fb=feval(f,b); it=0; whileabs(b-a)>tol, it=it+1; x=a/2+b/2; fx=feval(f,x); ifsign(fx)==sign(fa), a=x;fa=fx; else b=x;fb=fx; end end 11.牛顿法: 1.先定义函数后进行编程的方法 先需要定义不动点函数 需要计算的函数f functionv=f(x) v=x^5-3*x-10; end 需要计算的函数的导数g functionv=g(x) v=5*x^4-3; end 2.再定义编程: functionv=newton(a,f,g,maxit,tol) it=0; x=a; whileit<=maxit&abs(feval(f,x))>tol, it=it+1; x=x-feval(f,x)/feval(g,x); end v=[x,it]; end 12.牛顿下山法: 1.先定义函数后进行编程的方法 先需要定义不动点函数 需要计算的函数f functionv=f(x) v=x^2+sin(10*x)-1 end 需要计算的函数的导数g functionv=g(x) v=2*x+10*cos(10*x) end 2.再定义编程1: functionv=newtonxiashang(x0,f,g,maxit,tol) x=x0; it=0; whileit<=maxit&abs(feval(f,x))>tol, it=it+1; d=-feval(f,x)/feval(g,x); lambda=1; isdone=0; while~isdone, xn=x+lambda*d; ifabs(feval(f,xn)) isdone=1; elselambda=lambda*0.5; end end x=xn; end v=[x,it]; end 3.再定义编程2: functionv=newtonxiashang2(x0,f,g,maxit,tol) x=x0; it=0; whileit<=maxit&abs(feval(f,x))>tol, it=it+1; d=-feval(f,x)/feval(g,x); lambda=1; whileabs(feval(f,x+lambda*d))>=abs(feval(f,x)), lambda=0.5*lambda; end x=x+lambda*d; end v=[x,it]; end 13.割线法 1.先定义函数后进行编程的方法 先需要定义函数 需要计算的函数f functionv=f(x) v=x^5-3*x-10; end 2.再定义编程: functionv=gexian(a,b,f,maxit,tol) it=0; x0=a; x=b; whileit<=maxit&abs(feval(f,x))>tol, it=it+1; xt=x-feval(f,x)/(feval(f,x)-feval(f,x0))*(x-x0); x0=x; x=xt; end v=[x,it]; end 输入: >>gexian(1,2,'f',100,1e-5) ans= 1.72267.0000 14.乘幂法 function[t,y]=chenmifa(a,xinit,ep) v0=xinit; [tvti]=max(abs(v0)); lam0=v0(ti); u0=v0/lam0; flag=0; while(flag==0) v1=a*u0; [tvti]=max(abs(v1)); lam1=v1(ti); u0=v1/lam1; if(abs(lam0-lam1)<=ep) flag=1; end lam0=lam1; end t=lam1; y=u0; end 运行及其结果: >>a=[2,3,2;10,3,4;3,6,1]; >>xinit=[111]'; >>ep=0.0001; >>[ty]=chenmifa(a,xinit,ep) t= 11.0000 y= 0.5000 1.0000 0.7500 15.反幂法 function[ty]=fanmifa(a,xinit,ep) v0=xinit; [tvti]=max(abs(v0)); lam0=v0(ti); u0=v0/lam0; flag=0; while(flag==0) v1=a^-1*u0; [tvti]=max(abs(v1)); lam1=v1(ti); u0=v1/lam1; if(abs(lam0^-1-lam1^-1)<=ep) flag=1; end lam0=lam1; end t=lam1^-1; y=u0; end 运行提示: >>a=[126-6;6162;-6216]; >>xinit=[1-0.50.5]'; >>[ty]=fanmifa(a,xinit,0.0001) t= 4.4560 y= 1.0000 -0.6287 0.6287 16.结合原点平移反幂法 function[ty]=yuandian(a,p,xinit,ep) [nn]=size(a); v0=xinit; [tvti]=max(abs(v0)); lam0=v0(ti); u0=v0/lam0; flag=0; while(flag==0) v1=(a-p*eye(n))^-1*u0; [tvti]=max(abs(v1)); lam1=v1(ti); u0=v1/lam1; if(abs(lam0^-1-lam1^-1)<=ep) flag=1; end lam0=lam1; end t=lam1^-1+p; y=u0; end 运行结果: >>a=[621;231;111]; >>p=6; >>ep=0.0001; >>xinit=[111]'; >>[ty]=yuandian(a,p,xinit,ep) t= 7.2880 y= 1.0000 0.5229 0.2422 17.欧拉公式初值问题 1.先定义导数函数f: functionv=f(x,y) v=-y+x+1; end 2.定义欧拉公式编程: function[x,y]=oula(f,y0,a,b,n) y (1)=y0; h=(b-a)/n; x=a: h: b; fori=1: n; y(i+1)=y(i)+h*feval(f,x(i),y(i)); end end 18.改进的欧拉公式初值问题 1.先定义导数函数f: functionv=f(x,y) v=-y+x+1; end 要先运行一下。 2.定义欧拉公式编程: function[x,y]=gaijingoula(f,y0,a,b,n) y (1)=y0; h=(b-a)/n; x=a: h: b; fori=1: n; yp=y(i)+h*feval(f,x(i),y(i)); yc=y(i)+h*feval(f,x(i+1),yp); y(i+1)=0.5*(yc+yp); end end 19.梯形公式初值问题 1.先定义导数函数f: functionv=f(x,y) v=-y+x+1; end 要先运行一下。 2.定义欧拉公式编程: function[x,y]=tixing(f,y0,a,b,n) y (1)=y0; h=(b-a)/n; x=a: h: b; fori=1: n; y(i+1)=1/(0.5*h+1)*(y(i)+0.5*h*(feval(f,x(i),y(i))+x(i+1)+1)); end end 20.标准四阶四段龙格库塔公式初值问题 1.先定义导数函数f: functionv=f(x,y) v=y+x; end 要先运行一下。 2.定义欧拉公式编程: function[x,y]=longgekuta(f,y0,a,b,n) y (1)=y0; h=(b-a)/n; x=a: h: b; fori=1: n, k1=h*feval(f,x(i),y(i)); k2=h*feval(f,x(i)+0.5*h,y(i)+0.5*k1); k3=h*feval(f,x(i)+0.5*h,y(i)+0.5*k2); k4=h*feval(f,x(i)+h,y(i)+k3); y(i+1)=y(i)+(k1+2*k2+2*k3+k4)/6; end end 21.基本问题 1.特征值: eig(A) 2.数值显示位数: >>formatlong >>formatshort 3.左除和右除: AX=B,则X=A\B;左除\ XB=A,则X=A/B;右除/(右除是真除) 4.特殊符号 In(x)=log(x) 根号用sqrt(x) 5.绘图功能 Plot(x,y,style) >>x=0: 0.5: 10; >>y=x; >>holdon; >>plot(x,y,’g+: ’) 后面style的意义: y,g,r,b分别指颜色黄色、绿色、红色、蓝色; o*+p代表节点的形状,圆形、*号、+号、五角星; --、-、: 代表线型,分别为虚线、线段、点线
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 分析 编程 总结
![提示](https://static.bingdoc.com/images/bang_tan.gif)