zoukankan      html  css  js  c++  java
  • 「学习记录」《数值分析》第三章计算实习题(Python语言)

    第三题暂缺,之后补充。

    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')
    
    如非注明,原创内容遵循GFDLv1.3发布;其中的代码遵循GPLv3发布。
  • 相关阅读:
    div中嵌套div中使用margin-top失效问题
    thinkphp点击导航变色
    thinkphp I() 方法
    判断是手机端还是电脑端 isMobile()
    手机端H5 header定义样式
    AR.Drone 2.0四轴飞机体验:最好的玩具航拍器
    这是一个专注于电脑技术、软件应用、互联网、嵌入式,电子技术行业等的原创IT博客
    ul li列子
    [HTML]去除li前面的小黑点,和ul、LI部分属性
    Bad update sites
  • 原文地址:https://www.cnblogs.com/samhx/p/9652075.html
Copyright © 2011-2022 走看看