精华内容
下载资源
问答
  • 数值分析计算实习题
    2021-04-24 10:48:25

    《数值分析计算实习题》由会员分享,可在线阅读,更多相关《数值分析计算实习题(17页珍藏版)》请在人人文库网上搜索。

    1、数值分析计算实习题姓名: 学号: 班级: 第二章1、程序代码Clear;clc;x1=0.2 0.4 0.6 0.8 1.0;y1=0.98 0.92 0.81 0.64 0.38;n=length(y1);c=y1(:);for j=2:n %求差商for i=n:-1:jc(i)=(c(i)-c(i-1)/(x1(i)-x1(i-j+1);endendsyms x df d;df(1)=1;d(1)=y1(1);for i=2:n %求牛顿差值多项式df(i)=df(i-1)*(x-x1(i-1);d(i)=c(i-1)*df(i);endP4=vpa(sum(d),5) %P4即为4次牛。

    2、顿插值多项式,并保留小数点后5位数pp=csape(x1,y1, variational);%调用三次样条函数q=pp.coefs;q1=q(1,:)*(x-.2)3;(x-.2)2;(x-.2);1;q1=vpa(collect(q1),5)q2=q(1,:)*(x-.4)3;(x-.4)2;(x-.4);1;q2=vpa(collect(q2),5)q3=q(1,:)*(x-.6)3;(x-.6)2;(x-.6);1;q3=vpa(collect(q3),5)q4=q(1,:)*(x-.8)3;(x-.8)2;(x-.8);1;q4=vpa(collect(q4),5)%求解并化简多项式2。

    3、、运行结果P4 =0.98*x - 0.3*(x - 0.2)*(x - 0.4) - 0.625*(x - 0.2)*(x - 0.4)*(x - 0.6) - 0.20833*(x - 0.2)*(x - 0.4)*(x - 0.8)*(x - 0.6) + 0.784q1 =- 1.3393*x3 + 0.80357*x2 - 0.40714*x + 1.04q2 =- 1.3393*x3 + 1.6071*x2 - 0.88929*x + 1.1643q3 =- 1.3393*x3 + 2.4107*x2 - 1.6929*x + 1.4171q4 =- 1.3393*x3 + 3.2。

    4、143*x2 - 2.8179*x + 1.86293、问题结果4次牛顿差值多项式= 0.98*x - 0.3*(x - 0.2)*(x - 0.4) - 0.625*(x - 0.2)*(x - 0.4)*(x - 0.6) - 0.20833*(x - 0.2)*(x - 0.4)*(x - 0.8)*(x - 0.6) + 0.784三次样条差值多项式第三章1、程序代码Clear;clc;x=0 0.1 0.2 0.3 0.5 0.8 1;y=1 0.41 0.5 0.61 0.91 2.02 2.46;p1=polyfit(x,y,3)%三次多项式拟合p2=polyfit(x,y,4)。

    5、%四次多项式拟合y1=polyval(p1,x);y2=polyval(p2,x);%多项式求值plot(x,y,c-,x,y1,r:,x,y2,y-.)p3=polyfit(x,y,2)%观察图像,类似抛物线,故用二次多项式拟合。y3=polyval(p3,x);plot(x,y,c-,x,y1,r:,x,y2,y-.,x,y3,k-)%画出四种拟合曲线2、运行结果p1 =-6.6221 12.8147 -4.6591 0.9266p2 =2.8853 -12.3348 16.2747 -5.2987 0.9427p3 =3.1316 -1.2400 0.73563、问题结果三次多项式拟合P。

    6、1=四次多项式拟合P2=二次多项式拟合P3=第四章1、 程序代码1)建立函数文件f.m:function y=fun(x);y=sqrt(x)*log(x);2)编写程序:a. 利用复化梯形公式及复化辛普森公式求解:Clear;clc;h=0.001;%h为步长,可分别令h=1,0.1,0.01,0.001n=1/h;t=0;s1=0;s2=0;for i=1:n-1t=t+f(i*h);endT=h/2*(0+2*t+f(1);T=vpa(T,7) %梯形公式for i=0:n-1s1=s1+f(h/2+i*h);endfor i=1:n-1s2=s2+f(i*h);endS=h/6*(0+。

    7、4*s1+2*s2+f(1);S=vpa(S,7) %辛普森公式a复化梯形公式和复化辛普生公式程序代码也可以是:Clear;clc;x=0:0.001:1; %h为步长,可分别令h=1,0.1,0.01,0.001y=sqrt(x).*log(x+eps);T=trapz(x,y);T=vpa(T,7)(只是h=1的运行结果不一样,T=1.*10(-16),而其余情况结果都相同)Clear;clc;f=inline(sqrt(x).*log(x),x);S=quadl(f,0,1);S=vpa(S,7)b. 利用龙贝格公式求解:Clear;clc;m=14;%m+1即为二分次数h=2;for 。

    8、i=1:mh=h/2;n=1/h;t=0;for j=1:n-1t=t+f(j*h);endT(i)=h/2*(0+2*t+f(1);%梯形公式endfor i=1:m-1for j=m:i+1T(j)=4i/(4i-1)*T(j)-1/(4i-1)*T(j-1);%通过不断的迭代求得T(j),即T表的对角线上的元素。endendvpa(T(m),7)2、运行结果T =-0.S =-0.ans =-0.3、问题结果a. 利用复合梯形公式及复合辛普森公式求解:步长h10.10.010.001梯形求积T=01.*10(-16)-0.-0.-0.辛普森求积S=-0.-0.-0.-0.b. 利用龙贝格。

    9、公式求解:通过15次二分,得到结果:-0.第五章1、程序代码(1)LU分解解线性方程组:Clear;clc;A=10 -7 0 1-3 2. 6 25 -1 5 -12 1 0 2;b=8;5.;5;1;m,n=size(A); L=eye(n); U=zeros(n); flag=ok; for i=1:n U(1,i)=A(1,i); end for r=2:n L(r,1)=A(r,1)/U(1,1); end for i=2:n for j=i:n z=0; for r=1:i-1 z=z+L(i,r)*U(r,j); end U(i,j)=A(i,j)-z; end if abs(U。

    10、(i,i)k temp=Aug(k,:); Aug(k,:)=Aug(r,:); Aug(r,:)=temp; end if Aug(k,k)=0, error(对角元出现0), end % 把增广矩阵消元成为上三角for p = k+1:n Aug(p,:)=Aug(p,:)-Aug(k,:)*Aug(p,k)/Aug(k,k); end end % 解上三角方程组A = Aug(:,1:n); b = Aug(:,n+1); x(n) = b(n)/A(n,n); for k = n-1:-1:1 x(k)=b(k); for p=n:-1:k+1 x(k) = x(k)-A(k,p)*x。

    11、(p); end x(k)=x(k)/A(k,k); enddetA=det(A)2、 运行结果1) LU分解解线性方程组L =1.0e+006 *0.0000 0 0 0-0.0000 0.0000 0 00.0000 -2.5000 0.0000 00.0000 -2.4000 0.0000 0.0000U =1.0e+007 *0.0000 -0.0000 0 0.00000 -0.0000 0.0000 0.00000 0 1.5000 0.57500 0 0 0.0000x =-0.0000-1.00001.00001.0000detA =-762.00012)列主元消去法detA 。

    12、=762.0001ans =0.0000-1.00001.00001.00003、 问题结果1) LU分解解线性方程组L=U=x=(0.0000,-1.0000,1.0000,1.0000)TdetA=-762.0012)列主元消去法x=(0.0000,-1.0000,1.0000,1.0000)TdetA=762.001第六章1、程序代码(1)Jacobi迭代Clear;clc;n = 6; %也可取n=8,10H = hilb(n); b = H * ones(n, 1); e = 0.00001; for i = 1:n if H(i, i)=0 对角元为零,不能求解 return en。

    13、dendx = zeros(n, 1); k = 0; kend = 10000; r = 1; while ke x0 = x; for i = 1:n s = 0; for j = 1:i - 1 s = s + H(i, j) * x0(j); endfor j = i + 1:n s = s + H(i, j) * x0(j); endx(i) = b(i) / H(i, i) - s / H(i, i); endr = norm(x - x0, inf); k = k + 1; endif kkend 迭代不收敛,失败 else 求解成功x k end(2)SOR迭代1)程序代码fu。

    14、nction s = SOR(n, w); H = hilb(n); b = H*ones(n, 1); e = 0.00001; for i = 1:n if H(i,i)=0 对角线为零,不能求解return endendx = zeros(n, 1); k = 0;kend = 10000; r = 1; while ke x0 = x; for i = 1:n s = 0; for j = 1:i - 1s = s + H(i, j) * x(j); endfor j = i + 1:n s = s + H(i, j) * x0(j); endx(i) = (1 - w) * x0(i。

    15、) + w / H(i, i) * (b(i) - s);endr = norm(x - x0, inf); k = k + 1; endif kkend 迭代不收敛,失败 else 求解成功 xend2)从命令窗口中分别输入:SOR(6,1)SOR(8,1)SOR(10,1)SOR(6,1.5)SOR(8,1.5)SOR(10,1.5)2、运行结果Jacobi迭代:ans =迭代不收敛,失败SOR迭代:第七章1、程序代码(1)不动点迭代法1)建立函数文件:g.mfunction f=g(x)f(1)=20/(x2+2*x+10);2)建立函数文件:bdd.mfunction y, n = bdd(x, eps) if nargin=1eps=1.0e-8;elseif nargin=eps)&(n bdd(0)n =25ans =1.3688(2) 牛顿迭代x=021.66671.59211.38951.26671.13731.13731.1373k =8x1 =1.13733、问题结果(1)不动点迭代法x=1.3688 n=25 收敛太慢(2)牛顿迭代初值取0迭代次数k=8时,kx(k)1221.31.41.51.61.71.81。

    更多相关内容
  • 数值分析计算实习题答案 篇一数值分析 (第五版 )计算实习题第三章 t第二次作业 题一 x=-1:0.2:1;y=1./(1+25*x42; f1=polyfit(x,y,3) f=poly2sym(f1) y1=polyval(f1,x) x2=linspace(-1,1,10) y2=interp1(x,y,x2) plot...
  • 数值分析中一些 三次样条,龙格库塔,阿当姆斯等多种算法的实现
  • 数值分析(第五版)部分上机实习题报告和源代码(Python),题目为第二章(插值法)、第六章(解线性方程组的迭代法)、第七章(非线性方程与方程组的数值解法)、第九章( 常微分方程初值问题数值解法)
  • 北航数值分析上机编程第二,用带双步移位的QR分解法求解矩阵的全部特征值,并求相应的特征向量.
  • 数值分析计算实习第二 QR分解法求解矩阵特征值 代码亲测可用
  • 数值计算方法上机实习题 n 1 x 1 设 In dx 0 5 x 1 1 由递推公式 I n 5I n 1 从 , I0 =0.1823 出发计算 I 20 n 1 1 2 I 20 =0 I20 =10000 , 用 I n 1 I n 1 计算 I 0 5 5n (1) 由 I 0 计算 I 20 (2) 由 I 20 计算 I ...
  • 数值分析》为HUST计算机学院研究生阶段必修课,里面有配套的PPT,以及大量的例题、习题及参考答案。
  • 数值分析实验(参考答案).doc
  • 数值分析计算实习

    2018-12-27 11:51:57
    颜庆津的数值分析计算实习第一,无需修改直接调试运行即可,仅供参考!
  • 北航数值分析计算实习大作业1 幂法反幂法求特征值
  • 北航数值分析上机编程第三,分片二次差值,曲面拟合,Newton迭代法求解非线性方程组的解
  • 数值分析实习题

    2017-04-16 11:02:25
    天津大学计算机学院 数值分析 ,包括实验报告和代码
  • 第三暂缺,之后补充。 import matplotlib.pyplot as plt import numpy as np import scipy.optimize as so import sympy as sp x = sp.symbols('x') def calculate(expr_i, expr_j, expr_value,expr_omega): ...

    第三题暂缺,之后补充。

    import matplotlib.pyplot as plt
    import numpy as np
    import scipy.optimize as so
    import sympy as sp
    
    x = sp.symbols('x')
    
    def calculate(expr_i, expr_j, expr_value,expr_omega):
        ans=0
        for cnt,v in enumerate(expr_value):
            if isinstance(expr_i,(type(x),type(x*x))):
                actual_expr_i=expr_i.subs(x, v[0])
            elif expr_i==1: # which means 1 or 0
                actual_expr_i=1
            else:
                actual_expr_i=v[1]
            if isinstance(expr_j,(type(x),type(x*x))):
                actual_expr_j=expr_j.subs(x, v[0])
            else: # which means 1
                actual_expr_j=1
            if type(expr_omega)==type(1):
                ans=ans+expr_omega*actual_expr_i*actual_expr_j
            else:
                ans=ans+expr_omega[cnt]*actual_expr_i*actual_expr_j
    
        return ans
    
    def least_squares(_psi,_value,_omega):
        G=np.empty((len(_psi),len(_psi)))
        d=np.empty(len(_psi))
        for i,expr_i in enumerate(_psi):
            for j,expr_j in enumerate(_psi):
                G[i][j]=calculate(expr_i,expr_j,_value,_omega)
            d[i]=calculate(0,_psi[i],_value,_omega)
        a=np.linalg.solve(G, d.T) # Oh, I love solve()!
        ls_f_x=0
        for i,element in enumerate(a):
            ls_f_x=ls_f_x+element*_psi[i]
        return ls_f_x    
    
    psi=[1,x,x*x,x*x*x]
    #psi=[1,x]
    value=np.loadtxt('input.txt')
    omega=1
    #omega=[2,1,3,1,1]
    ls_f_x=least_squares(psi,value,omega)
    print(ls_f_x)
    f_x=1/(1+25*x*x)
    # p1=sp.plot(f_x,ls_f_x,(x,-1,1))
    
    """
    # using system functions
    def func(p,x):
        a0,a1,a2,a3 = p
        return a0+a1*x+a2*x*x+a3*x*x*x
    def err(p,x,y):
        return func(p,x)-y
    arg_f=(np.array([x[0] for x in value[:,:1]]),np.array([y[0] for y in value[:,1:2]]))
    coef_ls = so.leastsq(err, [1,1,1,1], args=arg_f)
    print(coef_ls)
    system_ls_f_x=0
    for i,element in enumerate(coef_ls[0]):
        system_ls_f_x=system_ls_f_x+element*psi[i]
    print(system_ls_f_x)
    p1=sp.plot(f_x,ls_f_x,system_ls_f_x,(x,-1,1),show=False)
    p1[1].line_color='r'
    p1[2].line_color='g'
    p1.show()
    """
    
    # problem 2:
    fig=plt.figure()
    
    value=np.array([[0,1],[0.1,0.41],[0.2,0.50],[0.3,0.61], [0.5,0.91],[0.8,2.02],[1.0,2.46]])
    
    ls_f_x=least_squares(psi,value,omega)
    print_f=sp.lambdify(x,ls_f_x,"numpy")
    print_x=np.linspace(-1,1,100)
    print_y=print_f(print_x)
    plt.plot(print_x,print_y,label='x^3')
    
    psi=[1,x,x**2,x**3,x**4]
    ls_f_x=least_squares(psi,value,omega)
    print_f=sp.lambdify(x,ls_f_x,"numpy")
    print_y=print_f(print_x)
    plt.plot(print_x,print_y,label='x^4')
    
    psi=[1,x]
    ls_f_x=least_squares(psi,value,omega)
    print_f=sp.lambdify(x,ls_f_x,"numpy")
    print_y=print_f(print_x)
    plt.plot(print_x,print_y,label='kx+b')
    
    plt.scatter(np.array([x[0] for x in value[:,:1]]),np.array([y[0] for y in value[:,1:2]]),marker='x',label='data')
    plt.legend(loc='best')
    plt.savefig('output.svg')
    
    展开全文
  • 利用Python完成了第一章的计算实习题——现在能搜到的基本上都是MATLAB版,或者是各种零碎的版本。这里将它完整的po出来,也方便自己以后查阅,提高自己的Python水平和数值分析的功底。 代码如下: (第一题使用的...

    在假期利用Python完成了《数值分析》第二章的计算实习题,主要实现了牛顿插值法和三次样条插值,给出了自己的实现与调用Python包的实现——现在能搜到的基本上都是MATLAB版,或者是各种零碎的版本。
    代码如下:
    (第一题使用的自己的程序,第二第三题使用的Python自带库)

    import math
    
    import matplotlib.pyplot as plt
    import numpy as np
    import pandas as pd
    from numpy.linalg import solve
    from scipy import interpolate
    from scipy.interpolate import lagrange
    
    plt.rc('figure',figsize=(20,15))
    
    print("Problem I:")
    
    given_x=[0.2,0.4,0.6,0.8,1.0]
    given_y=[0.98,0.92,0.81,0.64,0.38]
    given_times=4
    x_range=(0,1.1,0.02)
    
    #@brief: Convert(begin,end,interval) to a list, but interval can be float numbers.
    def process_xpara(xpara):
        max_times=0
        if 0<xpara[2]<1:
            tmp_xpara_interval=xpara[2]
            while tmp_xpara_interval-int(tmp_xpara_interval)!=0:
                max_times=max_times+1
                tmp_xpara_interval=tmp_xpara_interval*10
        max_times=10**max_times
        return [i/max_times for i in range(int(xpara[0]*max_times),int(xpara[1]*max_times),int(xpara[2]*max_times))]
    
    
    def divide_difference(x,y,times):
        now=[(x[i],y[i]) for i in range(len(x))]
        ans=[now[0]]
        for order in range(1,times+1):
            tmp=[]
            for i in range(1,len(now)):
                tmp.append((x[order+i-1]-x[i-1],(now[i][1]-now[i-1][1])/(x[order+i-1]-x[i-1])))
            now=tmp
            ans.append(now[0])
        return ans
    
    def get_func_value_newton(xcoef,x,xorigin):
        ans=0
        for i in range(len(xcoef)):
            tmp=xcoef[i][1]
            for j in range(i):
                tmp=tmp*(x-xorigin[j])
            ans=ans+tmp
        return ans
    """
    #@param: xpara(xbegin,xend,xinterval) fpara[f[x_1~i]]
    """
    # spec_i=[0.2+0.08*x for x in (0,1,10,11)]
    
    def newton_interpolate(xpara,fpara,xorigin):
        x_discrete_value=process_xpara(xpara)
        return [(x,get_func_value_newton(fpara,x,xorigin)) for x in x_discrete_value]
    
    parameters=divide_difference(given_x,given_y,given_times)
    newton_interpolate_value=newton_interpolate(x_range,parameters,given_x)
    
    fig=plt.figure()
    sub_fig1=fig.add_subplot(2,2,1)
    sub_fig1.set_title("Problem I")
    sub_fig1.plot([var[0] for var in newton_interpolate_value],[var[1] for var in newton_interpolate_value],label='Newton')
    
    # l_f=lagrange(given_x,given_y)
    # tmpara=process_xpara(x_range)
    # plt.plot(tmpara,[l_f(x) for x in tmpara])
    
    # 三次样条插值
    n=len(given_x)
    h=[]
    f0p=0
    fnp=0
    for i in range(1,len(given_x)):
        h.append(given_x[i]-given_x[i-1])
    miu=[0] # 0 should not be used
    lam=[1]
    d=[6/h[0]*((given_y[1]-given_y[0])/(given_x[1]-given_x[0])-f0p)]
    for j in range(1,len(h)):
        miu.append(h[j-1]/(h[j-1]+h[j]))
        lam.append(h[j]/(h[j-1]+h[j]))
        d.append(6*((given_y[j+1]-given_y[j])/(given_x[j+1]-given_x[j])-(given_y[j-1]-given_y[j])/(given_x[j-1]-given_x[j]))/(h[j-1]+h[j]))
    miu.append(1)
    d.append(6/h[-1]*(fnp-(given_y[-1]-given_y[-2])/(given_x[-1]-given_x[-2])))
    
    A=np.zeros((n,n))
    for i in range(n):
        A[i][i]=2
        if i!=n-1:
            A[i][i+1]=lam[i]
        if i!=0:
            A[i][i-1]=miu[i]
    C=solve(A,np.array(d).T)
    # print(C)
    
    def get_func_value_cubic_spline(mtuple, xtuple, ytuple, x):
        return mtuple[0]/(6*(xtuple[1]-xtuple[0]))*(xtuple[1]-x)**3+mtuple[1]/(6*(xtuple[1]-xtuple[0]))*(x-xtuple[0])**3+(ytuple[0]-(mtuple[0]*(xtuple[1]-xtuple[0])**2/6))*(xtuple[1]-x)/(xtuple[1]-xtuple[0])+(ytuple[1]-(mtuple[1]*(xtuple[1]-xtuple[0])**2/6))*(x-xtuple[0])/(xtuple[1]-xtuple[0])
    
    def cubic_spline_interpolate(xpara, mpara, x, y):
        fun_value=[]
        x_discrete_value=process_xpara(xpara)
        for j in range(len(x)-1):
            ok_value=[(element,get_func_value_cubic_spline((mpara[j],mpara[j+1]),(x[j],x[j+1]),(y[j],y[j+1]),element)) for element in x_discrete_value if x[j]<=element<x[j+1]]
            fun_value=fun_value+ok_value
        return fun_value
    cubic_spline_interpolate_value=cubic_spline_interpolate(x_range,C.tolist(),given_x,given_y)
    
    sub_fig1.plot([var[0] for var in cubic_spline_interpolate_value],[var[1] for var in cubic_spline_interpolate_value],label='Cubic')
    
    sub_fig1.legend(loc='best')
    
    def get_func_x(x):
        return 1/(1+25*x*x)
    
    given_x=np.linspace(-1,1,10)
    given_y=get_func_x(given_x) #p.array([get_func_x(x) for x in given_x])
    display_x=np.linspace(-1,1,100)
    display_y=get_func_x(display_x)
    
    sub_fig2=fig.add_subplot(2,2,2)
    sub_fig2.set_title("Problem II(Alpha): Using System Functions")
    
    c_x=interpolate.interp1d(given_x, given_y, kind = 'cubic')
    l_x=lagrange(given_x,given_y)
    sub_fig2.plot(display_x,l_x(display_x))
    sub_fig2.plot(display_x,c_x(display_x))
    sub_fig2.plot(display_x,display_y)
    
    sub_fig3=fig.add_subplot(2,2,3)
    sub_fig3.set_title("Problem II(Beta): Using System Functions")
    
    given_x=np.linspace(-1,1,20)
    given_y=get_func_x(given_x) #p.array([get_func_x(x) for x in given_x])
    c_x=interpolate.interp1d(given_x, given_y, kind = 'cubic')
    l_x=lagrange(given_x,given_y)
    sub_fig3.plot(display_x,l_x(display_x))
    sub_fig3.plot(display_x,c_x(display_x))
    sub_fig3.plot(display_x,display_y)
    
    fig_problem_three=plt.figure()
    
    given_x=[0,1,4,9,16,25,36,49,64]
    given_y=[0,1,2,3,4,5,6,7,8]
    display_big_x=np.linspace(0,64,200)
    display_small_x=np.linspace(0,1,50)
    sub_fig4=fig_problem_three.add_subplot(2,1,1)
    l_x=lagrange(given_x,given_y)
    c_x=interpolate.interp1d(given_x, given_y, kind = 'cubic')
    sub_fig4.plot(display_big_x,l_x(display_big_x),label='Lagrange')
    sub_fig4.plot(display_big_x,c_x(display_big_x),label='Cubic')
    sub_fig4.plot(display_big_x,np.sqrt(display_big_x),label='Origin')
    sub_fig4.legend(loc='best')
    
    sub_fig5=fig_problem_three.add_subplot(2,1,2)
    sub_fig5.plot(display_small_x,l_x(display_small_x),label='Lagrange')
    sub_fig5.plot(display_small_x,c_x(display_small_x),label='Cubic')
    sub_fig5.plot(display_small_x,np.sqrt(display_small_x),label='Origin')
    sub_fig5.legend(loc='best')
    展开全文
  • 北京航空航天大学数值分析计算实习题目,所有习题,C#,VS2008环境,保证能用且答案正确
  • 西南交通大学数值分析上机实习题(内含各种问题源代码)
  • 数值分析李庆杨第五章习题第十八题计算实习题第三题分析结果计算实习题第四题第一问第二问收获 第十八题 学习条件数 逆矩阵 矩阵转置 特征向量、特征值的求法 X = [100,99;99,98]; Y = inv(X); B = transpose(X); A...

    第十八题

    学习条件数
    逆矩阵 矩阵转置
    特征向量、特征值的求法

    X = [100,99;99,98];
    Y = inv(X);
    B = transpose(X);
    A = B*A;
    %A=[19801,19602;19602,19405]
    e = eig(A);
    sqrt(e(2)/e(1));
    cond(X,2)
    cond(X,inf)
    

    计算实习题第三题

    题目:从理论结果和实际计算两方面分析线性方程组AX=B解的相对误差及相对误差的关系

    A = [10,7,8,7;7,5,6,5;8,6,10,9;7,5,9,10];
    b = [32;23;33;31];
    x0 = [1;1;1;1];
    A_det = det(A)
    A_eig = eig(A)
    A_cond2 = cond(A,2)
    A_cond1 = cond(A,1)
    delta_A = [0,0,0.1,0.2;0.08,0.04,0,0;0,-0.02,-0.11,0;-0.01,0,0,-0.02];
    x1 = (A+delta_A)\b;
    delta_x = x1-x0
    delta_xnorm2 = norm(delta_x,2)
    wucha_x = norm(delta_x,2)/norm(x0,2)
    wucha_A = norm(delta_A,2)/norm(A,2)
    A_ni = inv(A);
    condition1 = norm(A_ni,1)*norm(delta_A,1)
    condition2 = norm(A_ni*delta_A,1)
    
    

    A_det = 1.0000
    A_eig = 4×1
    0.0102
    0.8431
    3.8581
    30.2887

    A_cond2 = 2.9841e+03
    A_cond1 = 4.4880e+03
    delta_x = 4×1
    -10.5863
    17.3741
    -4.2258
    2.5240

    delta_xnorm2 = 20.9322
    wucha_x = 10.4661
    wucha_A = 0.0076
    condition1 = 29.9200
    condition2 = 16.8200

    分析结果

    从计算结果看,A的2范数误差比是0.0076,表示矩阵A的扰动很小。x结果的2范数误差比是10.4661,相比较与A的扰动,x结果的改变很大。这表示矩阵A对于误差相当敏感,A的微小改变就能引起解x较大的改变。结合理论计算,A的条件数cond(A)1的值为4.4880e+03,远远大于1,它是病态的,如果一个系统是病态 的,就很难用一般的计算方法。联系书本上168面推导的相对误差的上界公式,可以明显看到,对于病态方程,||A-1*δA|| = 29.92>1,另一条件||A-1|| ||δA|| = 16.82>1,无法使用公式求出相对误差的上界。

    计算实习题第四题

    题目:希尔伯特矩阵

    第一问

    a = zeros(6,1);
    for i = 2:size(a)
    M = ones(i, 1)*(1:i);
    N = M';
    H = 1 ./ (M+N-1);
    a(i) = cond(H,inf);
    end
    axis equal
    plot(1:6,a,'r-','LineWidth',1.5);
    xlabel('迭代次数');
    ylabel('无穷条件数');
    
    迭代次数 ncond(Hn)
    11
    227
    3748
    42.8375e+04
    59.4366e+05
    62.9070e+07

    第二问

    format long
    a = zeros(13,1);
    for i = 1:13
        M = ones(i, 1)*(1:i);
        N = M';
        H = 1 ./ (M+N-1);
        x = ones(i,1);
        x_new = GaussElimination(H,(H*x)') 
        delta_x = x-x_new
        rn = x - H*x_new
        chazhi = max(abs(delta_x));
        a(i) = cond(H,inf);
    end
    plot(1:13,a,'r-','LineWidth',1.5);
    xlabel('迭代次数');
    ylabel('无穷条件数');
    
    迭代次数 n有效位数max(Δx)
    2156.6613e-16
    3141.2546e-14
    4125.9797e-13
    5122.5913e-12
    6104.2542e-10
    781.5222e-08
    874.8704e-07
    951.6869e-05
    1044.3691e-04
    1120.0097
    1210.2042
    1306.4935

    收获

    script 的使用
    不要重复运行

    还不会matlab输出 %d i 的那种表达

    使用福昕pdf编辑器
    手动 可以动图
    文本可以写字
    右键插入文件 即可拼接

    展开全文
  • 数值分析大作业 第一 学号 姓名 学院 授课教师 2016.10.24 一算法设计方案 1求1501 a)利用幂法计算出矩阵A按模最大的特征值设其为m b)令矩阵B=A-mI c)令m'=m'+m由计算过程可知 d) 利用反幂法计算s其中通过使用...
  • 研究生数值计算实习大作业第一次和第二次,作业文档,内有代码,可编译运行
  • 数值分析思考9.doc

    2020-02-15 01:46:58
    数值分析思考9 一个算法局部误差和整体误差的区别是什么如何定义常微分方程数值方法的阶 称 为某方法在点的整体截断误差设是准确的用某种方法计算时产生的截断误差称为该方法的局部截断误差可以知道整体误差来自于...
  • 2010年数值分析大作业第二,希望在学习数值分析上能够帮学弟学妹一些忙,更好的进步~
  • 北京航空航天大学研究生数学必修课程,数值积分大作业。
  • 插值法 下列数据点的插值 x 0 1 4 9 16 25 36 49 64 y 0 1 2 3 4 5 6 7 8 可以得到平方根函数的近似在区间[0,64]上作图. (1)用这9个点作8次多项式插值Ls(x. (2)用三次样条(第一边界条件)程序求S(x....
  • 数值分析第四章上机实习题

    千次阅读 2013-10-31 21:02:52
    %考虑到会有NaN的情况,这句是把所有NaN的位置找出来,例如第一中左端点值 for ix = 1:length(findnan)%然后用极限进行替换,第一可以直接换成0. v(findnan(ix)) = limit(fx,u(findnan(ix))); end I = (v(1) + v...
  • 数值分析第九章作业

    2020-12-30 08:35:47
    第九章 欧拉法 function [X,Y] = Euler(a,b,h,f) n = abs((a-b)/h); X =a:h:b; Y =zeros(1,n); Y(1) = 1; for i = 1:n Y(i+1) = Y(i) +h* f(X(i),Y(i)); end 改进的欧拉法 function [X,Y] = Euler_gaijin...for i =
  • 数值分析第五版答案,第一章 绪论 第二章 插值法 第三章 函数逼近与曲线拟合
  • 北航数值分析上机编程第一,求矩阵的特征值、条件数、行列式,用到幂法、反幂法,还涉及求解线性方程组等。
  • 数值分析》所有算法、例题、课后的MATLAB程序实现 并且每个程序都附有详尽的说明和例子,绝对吐血珍藏 《数值分析》所有算法、例题、课后的MATLAB程序实现 并且每个程序都附有详尽的说明和例子,绝对吐血珍藏
  • 数值分析(数值计算) 数学建模 实验报告 matlab程序

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 2,114
精华内容 845
关键字:

数值分析计算实习题

友情链接: ariglxr.rar