• MATLAB数值分析实验


    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

     

  • 相关阅读:
    01 mybatis框架整体概况(2018.7.10)-
    第一课(2018.7.10)
    JavaEE 企业级分布式高级架构师课程_汇总贴
    5-1条件运算符 & 5-2
    5-3运算符的优先级
    4-3逻辑非运算符及案例 & 4-4
    4-1逻辑与运算符介绍 & 4-2逻辑或运算符介绍
    3-3if-else条件结构 & 3-4 & 3-5
    3-2if条件结构
    3-1关系运算符
  • 原文地址:https://www.cnblogs.com/kailugaji/p/6973002.html
Copyright © 2020-2023  润新知