第三题暂缺，之后补充。

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完成了《数值分析》第二章的计算实习题，主要实现了牛顿插值法和三次样条插值，给出了自己的实现与调用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')

第十八题

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

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)₁的值为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('无穷条件数');

第二问

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('无穷条件数');

收获

script 的使用
不要重复运行

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

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