1.用Newton迭代法求方程 的第一个正根.
作者:凯鲁嘎吉 - 博客园http://www.cnblogs.com/kailugaji/
newton.m: function x1=newton(x0,eps) format long format compact x1=x0-dao(x0); while abs(x1-x0)>eps x0=x1; x1=x0-dao(x0); end dao.m: function y=dao(x) y=tan(x)-exp(x); y1=tan(x)^2 - exp(x) + 1; y=y/y1;
结果:
>> x1=newton(1,1e-6)
x1 =
1.306326940423080
2.作矩阵 的LU分解.
lu12.m: function [l,u]=lu12(a,n) for k=1:n-1 for i=k+1:n a(i,k)=a(i,k)/a(k,k); for j=k+1:n a(i,j)=a(i,j)-a(i,k)*a(k,j); end end end l=eye(n); u=zeros(n,n); for k=1:n for i=k:n u(k,i)=a(k,i); end end for k=1:n for j=1:k-1 l(k,j)=a(k,j); end end 结果: >> a=[4 1 1 1;8 5 1 3;12 -3 7 2;4 10 2 7]; >> [l,u]=lu12(a,4) l = 1 0 0 0 2 1 0 0 3 -2 1 0 1 3 2 1 u = 4 1 1 1 0 3 -1 1 0 0 2 1 0 0 0 1
3.用Jacobi迭代法求解方程组, 其中.
jacobi.m: function x=jacobi(a,b,x0,n,tol,m) x=zeros(n,1); for k=0:m for i=1:n s=0; for j=1:n if j~=i s=s+a(i,j)*x0(j,1); end end x(i,1)=(b(i,1)-s)/a(i,i); if norm(x-x0,inf)<tol break; end x0(i,1)=x(i,1); end end 结果: >> a=[4 1 -1 1;-1 4 -1 1;1 2 5 -1;3 2 -1 7]; >> b=[2 8 8 10]'; >> x0=[0 0 0 0]'; >> x=jacobi(a,b,x0,4,1e-6,50) x = -0.000000983453000 2.000001278748222 0.999997309599650 0.999999964663427
4.用复化的辛甫生方法计算.
simpson.m: function [SI,Y,esp]=simpson(a,b,m) %a,b为区间左右端点,xps(x)为求积公式,m*2等分区间长度 h=(b-a)/(2*m); SI0=xps(a)+xps(b); SI1=0; SI2=0; for i=1:((2*m)-1) x=a+i*h; if mod(i,2)==0 SI2=SI2+xps(x); else SI1=SI1+xps(x); end end SI=vpa(h*(SI0+4*SI1+2*SI2)/3,10); syms x Y=vpa(int(xps(x),x,a,b),10); esp=abs(Y-SI); xps.m: function y=xps(x) y=exp(x^2)-sin(x)/x; 结果: >> [SI,Y,esp]=simpson(1,3,10) SI = 1443.251264 Y = 1442.179902 esp = 1.0713621257845176160117262043059
5.用改进的尤拉法解方程
euler22.m: function [B1,B2]=euler22(a,b,n,y0) %欧拉法解一阶常微分方程 %初始条件y0 h = (b-a)/n; %步长h %区域的左边界a %区域的右边界b x = a:h:b; m=length(x); %改进欧拉法 y = y0; for i=2:m y(i)=y(i-1)+h/2*( oula2(x(i-1),y(i-1))+oula2(x(i),y(i-1))+h*(oula2(x(i-1),x(i-1)))); B1(i)=x(i); B2(i)=y(i); end plot(x,y,'m-'); hold on; %精确解用作图 xx = x; f = dsolve('Dy=exp(x-y)+(x^2)*exp(-y)','y(0)=0','x');%求出解析解 y = subs(f,xx); %将xx代入解析解,得到解析解对应的数值 plot(xx,y,'k--'); legend('改进欧拉法','解析解'); oula2.m: function f=oula2(x,y) f=exp(x-y)+(x^2)*exp(-y); 结果: >> [B1,B2]=euler22(0,1,10,0) B1 = Columns 1 through 7 0 0.100000000000000 0.200000000000000 0.300000000000000 0.400000000000000 0.500000000000000 0.600000000000000 Columns 8 through 11 0.700000000000000 0.800000000000000 0.900000000000000 1.000000000000000 B2 = Columns 1 through 7 0 0.110758545903782 0.222173861791736 0.335492896789537 0.451351722029268 0.569931474513367 0.691088488902808 Columns 8 through 11 0.814464555075657 0.939577860819895 1.065894210026593 1.192879090561291
6.(1) 用拟合下列数据:
x |
2.36 |
3.73 |
5.951 |
8.283 |
f(x) |
14.1 |
16.2 |
18.3 |
21.4 |
LSM1.m: function [a,b,c]=LSM1(x,y,m) %x,y为序列长度相等的数据向量,m为拟合多项式次数 format short; A=zeros(m+1,m+1); for i=0:m for j=0:m A(i+1,j+1)=sum(x.^(i+j)); end b(i+1)=sum(x.^i.*y); end a=A'; p=fliplr(a'); %y=p[0]*x^m+p[1]*x^(m-1)+...+p[m-1]*x+p[m]; a=p(3); b=p(2); c=p(1); 结果: >> x=[2.36 3.73 5.951 8.283]; >> y=[14.1 16.2 18.3 21.4]; >> [a,b,c]=LSM1(x,y,2) a = 11.4457 b = 1.1866 c = 8.1204e-04
(2) 按如下插值原则,求Newton插值多项式:
x |
2.36 |
3.73 |
5.951 |
8.283 |
f(x) |
14.1 |
16.2 |
18.3 |
21.4 |
说明:最后,一定给清楚各多项式的系数!
newploy.m: function [A,C,L,wcgs,Cw]= newploy(X,Y) n=length(X); A=zeros(n,n); A(:,1)=Y'; q=1.0; c1=1.0; for j=2:n for i=j:n A(i,j)=(A(i,j-1)- A(i-1,j-1))/(X(i)-X(i-j+1)); end b=poly(X(j-1));q1=conv(q,b); c1=c1*j; q=q1; end C=A(n,n); b=poly(X(n)); q1=conv(q1,b); for k=(n-1):-1:1 C=conv(C,poly(X(k))); d=length(C); C(d)=C(d)+A(k,k); end L(k,:)=poly2sym(C); Q=poly2sym(q1); syms M wcgs=M*Q/c1; Cw=q1/c1; 结果: >> x=[2.36 3.73 5.951 8.283]; >> y=[14.1 16.2 18.3 21.4]; >> [A,C,L,wcgs,Cw]= newploy(x,y) A = 14.1000 0 0 0 16.2000 1.5328 0 0 18.3000 0.9455 -0.1636 0 21.4000 1.3293 0.0843 0.0418 C = 0.0418 -0.6674 4.4138 6.8506 L = (3015319848353441*x^3)/72057594037927936 - (3005803726105311*x^2)/4503599627370496 + (4969523982821561*x)/1125899906842624 + 7713109820116169/1125899906842624 wcgs = (M*(x^4 - (5081*x^3)/250 + (1273498286182623*x^2)/8796093022208 - (7485266609524121*x)/17592186044416 + 7633404131354389/17592186044416))/24 Cw = 0.0417 -0.8468 6.0325 -17.7287 18.0795 newpoly2.m: function y= newpoly2(X,Y,x) n=length(X); m=length(x); for t=1:m z=x(t); A=zeros(n,n);A(:,1)=Y'; q1=1.0; c1=1.0; for j=2:n for i=j:n A(i,j)=(A(i,j-1)- A(i-1,j-1))/(X(i)-X(i-j+1)); end q1=abs(q1*(z-X(j-1)));c1=c1*j; end C=A(n,n);q1=abs(q1*(z-X(n))); for k=(n-1):-1:1 C=conv(C,poly(X(k)));d=length(C); C(d)=C(d)+A(k,k); end y(k)= polyval(C, z); end 结果: >> y= newpoly2(x,y,15) y = 64.1181