小白学计算方法/数值计算 实验二

计算方法/数值计算 实验二

问题重述

问题分析

对于最高次数为n的拟合多项式，有：
$$f(x)=\sum _{i=0}^n{a}_i^{\ast }{h}_i(x)$$
其中：
$$h_i(x)=(x^{\ast }{)}^i$$
相应的法方程组为：

求解对应的参数向量即可。
此题对应的经验公式形式已给出，相应的法方程组形式也由此确定。本题中n=1。

实现代码

MATLAB实现

MATLAB代码

clc;
clear;

x_data = [1 2 4 8 16 32 64];
y_data = [4.22 4.02 3.85 3.59 3.44 3.02 2.59];

trans_x = log(x_data);
trans_y = log(y_data);

equations_matrix_x = [length(trans_x) sum(trans_x);
                      sum(trans_x) sum(trans_x .* trans_x)];
equations_matrix_y = [sum(trans_y);
                      sum(trans_x .* trans_y)];
                  
solved_parameters = equations_matrix_x ^ -1 * equations_matrix_y;

c = exp(solved_parameters(1));
lambda = solved_parameters(2);

disp("c: " + c)
disp("lambda: " + lambda)

plot(x_data, y_data, 'r*')
hold on
syms x
fplot(c*x.^lambda, 'b')

MATLAB实现结果

$$\begin{gathered} c:4.394 \\ lambda:-0.11074 \end{gathered}$$

Python实现

Python代码1

1. 计算主文件

# -*- coding: utf-8 -*-
"""
Created on Tue Oct 13 12:09:27 2020

@author: yfcod
"""

import numpy as np
import pandas as pd
import math


from solve_equations import sympy_solver

x_data = [1, 2, 4, 8, 16, 32, 64]
y_data = [4.22, 4.02, 3.85, 3.59, 3.44, 3.02, 2.59]

def trans_func(x):
    return math.log(x)

def reverse_trans_func(x):
    return math.pow(math.e, x)

def test_func(ln_x, test_c, test_lambda):
    ln_omega = test_lambda * ln_x + test_c
    return ln_omega

def trans_variable(data):
    trans_data = []
    for i in range(len(data)):
        trans_data.append(trans_func(data[i]))
    return trans_data

def generate_equations_matrix_x(x_data):
    return np.matrix([[len(x_data), np.sum(x_data)],
                      [np.sum(x_data), np.sum([x**2 for x in x_data])]])

def generate_equations_matrix_y(x_data, y_data):
    return np.array([np.sum(y_data),
                     np.sum([x_data[i] * y_data[i] for i in range(len(x_data))])])

trans_x = [trans_func(x) for x in x_data]
trans_y = [trans_func(y) for y in y_data]

equations_matrix_x = generate_equations_matrix_x(trans_x)
equations_matrix_y = generate_equations_matrix_y(trans_x, trans_y)

equations = [f'{equations_matrix_x[0, 0]} * a + {equations_matrix_x[0, 1]} * b - {equations_matrix_y[0]}',
             f'{equations_matrix_x[1, 0]} * a + {equations_matrix_x[1, 1]} * b - {equations_matrix_y[1]}']

parameters = ['a', 'b']

new_solver = sympy_solver(','.join(equations), ','.join(parameters))

new_solver.print_equations()
new_solver.print_parameters()
solved_parameters = list(new_solver.solve_basic_equations().values())

final_parameters = dict(zip(['c', 'lambda'], 
                            [reverse_trans_func(solved_parameters[0]), float(solved_parameters[1])]))

print('The final parameters is: ')
print(pd.DataFrame(final_parameters, index=[0]))

2. 基于sympy的方程求解工具类

# -*- coding: utf-8 -*-
"""
Created on Thu Jul 30 17:12:51 2020

@author: yfcod
"""

import sympy as sp


class sympy_solver(object):
    def __init__(self, equations, parameters):
        self.equations = equations
        self.parameters = parameters
        
        self.equations_list = self.string2list(self.equations)
        self.parameters_list = self.string2list(self.parameters)
        
        
    def __str__(self):
        return 'Please input the parameters split by \',\', input the equations in currect form'
    
    def print_equations(self):
        print('The equations:')
        print([item + ' = 0' for item in self.equations_list])
        
    def print_parameters(self):
        print('The parameters:')
        print(self.parameters_list)
    
    
    def string2list(self, string):
        if not isinstance(string, str):
            print('The parameters are not a string, input them again.')
            return
        
        new_list = string.split(',')
        
        for item_index, item in enumerate(new_list):
            new_list[item_index] = item.strip()
        
        return new_list
        
    
    def solve_basic_equations(self):
        para_symbol_list = []
        
        for item in self.parameters_list:
            para_symbol_list.append(sp.Symbol(item))
        
        solved_values = sp.solve(self.equations_list, para_symbol_list)
        
        print('solved values:')
        print(solved_values)
        
        print('format print:')
        sp.pprint(solved_values, use_unicode=True)
        
        return solved_values

这是一开始写的代码，对于法方程组的求解的处理不够精简（这里还用到了本人之前闲来无事用作消遣写的一个求解方程的工具类）。对于法方程组的求解完全可以用矩阵运算解决。并且该代码可重用性低，封装性差，效率不高。

Python代码1实现结果

The equations:
['7.0 * a + 14.556090791758852 * b - 8.74972855996226 = 0', '14.556090791758852 * a + 43.72122426655632 * b - 16.70484850951222 = 0']
The parameters:
['a', 'b']
solved values:
{a: 1.48023089175305, b: -0.110736303130349}
format print:
{a: 1.48023089175305, b: -0.110736303130349}
The final parameters is: 
         c    lambda
0  4.39396 -0.110736

Python代码2

为了封装成解决最小二乘法拟合多项式的工具类，增强代码的可重用性，故对代码1进行修改。

# -*- coding: utf-8 -*-
"""
Created on Tue Oct 13 14:50:40 2020

@author: yfcod
"""

import numpy as np
import pandas as pd
import math


class LeastSquareMethod():
    def __init__(self, x_data, y_data, max_dim):
        if len(x_data) != len(y_data):
            raise Exception('x_data\'s length can\'t match y_data\'s length')
        
        self.x_data = x_data
        self.y_data = y_data
        self.dim = max_dim
        
        self.equations_matrix_x = self.generate_matrix_x()
        self.equations_matrix_y = self.generate_matrix_y()
        
        self.parameters_list = self.solve_parameters()
        
    def generate_matrix_x(self):
        matrix_x = np.matrix([[np.sum([x**col for x in self.x_data]) for col in range(row, self.dim + row + 1)]
                             for row in range(self.dim + 1)
                             ])
        matrix_x[0, 0] = len(self.x_data)
        return matrix_x
    
    def generate_matrix_y(self):
        matrix_y = np.array([np.sum([(self.x_data[i]**row) * self.y_data[i] for i in range(len(self.x_data))])
                             for row in range(self.dim + 1)
                             ])
        return matrix_y
 
    def solve_parameters(self):
        parameters_matrix = np.matmul(self.equations_matrix_x.I, self.equations_matrix_y)
        parameters_list = parameters_matrix.tolist()[0]
        return parameters_list
 
def trans_func(x):
    return math.log(x)

def reverse_trans_func(x):
    return math.pow(math.e, x)    

x_data = [1, 2, 4, 8, 16, 32, 64]
y_data = [4.22, 4.02, 3.85, 3.59, 3.44, 3.02, 2.59]
max_dim = 1
    
trans_x = [trans_func(x) for x in x_data]
trans_y = [trans_func(y) for y in y_data]

test = LeastSquareMethod(trans_x, trans_y, max_dim)

final_parameters = dict(zip(['c', 'lambda'], 
                            [reverse_trans_func(test.parameters_list[0]), test.parameters_list[1]]))

print('The final parameters is: ')
print(pd.DataFrame(final_parameters, index=[0]))

我们构造LeastSquareMethod这个类，对其输入x样本数据，y样本数据，待拟合多项式最高次数。通过generate_matrix_x函数生成法方程组左边的系数矩阵，通过generate_matrix_y函数生成右边的向量，通过solve_parameters求解拟合多项式里的参数向量$[a_0,a_1,...,a_i]$。

Python代码2实现结果

The final parameters is: 
         c    lambda
0  4.39396 -0.110736

补充

对于曲线拟合问题，最高效率的工具当属MATLAB工具箱里的Curve Fitting Tool！