非线性方程求参数

发表于 | 更新于: | 分类于

Intro

最近遇到了一个问题，要求一个由多个非线性方程组成的方程的 参数
已经目标方程的形式 求解对应非线性方程的未知参数 例如：

$$ s(t)=a\cdot f1(t)+b\cdot f2(t)+c \cdot f(3t) $$

其中:

$$f1(t)= w \cdot sin(k\cdot t)$$

$$f2(t)= 1- \frac{1}{\sqrt1-w^2}\cdot e^{-\alpha}\cdot sin[(1-\sqrt{1-w^2})+arccos(w)]$$

$f(3)$是一个cdf 概率累计分布函数

Python 解法

首先找了个LMFIT的库,可以读文档使用

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
from numpy import sqrt, pi, exp, sin, arccos, cos
from lmfit import Model
from scipy import stats
import numpy as np
import matplotlib.pyplot as plt

__author__ = 'wangj'
__date__ = '2017/12/26 00:41'


def my_sin(x, w1, k):
    return w1 * sin(k * x)


def my_damping(x, w2, a, b, c):
    return w2 * (1 - 1 / sqrt(1 - a ** 2) * exp(-c) * sin(sqrt(1 - a ** 2) * b * x + arccos(a)))


def my_peak(x, m, n):
    return m * stats.poisson.cdf(x, n)


if __name__ == '__main__':
    x = np.linspace(1, 200, 100)
    y = 10 * cos(x) + sin(x)
    mod = Model(my_sin) + Model(my_damping)+ Model(my_peak)
    pars = mod.make_params(w1=0.0001, k=-0.4, w2=-0.001, a=0.1, b=1, c=-9,m=0,n=10)
    pars['a'].set(min=0.01, max=0.9999)
    result = mod.fit(y, pars, x=x)

    print(result.fit_report())

    plt.plot(x, y, 'bo')
    plt.plot(x, result.init_fit, 'k--')
    plt.plot(x, result.best_fit, 'r-')
    plt.show()

关键的问题在于初始params的选择，这边可以使用一下Matlab的lsqcurvefit函数，参考下初始值的设定

Matlab

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
t=1:15:400;
Ca = 10*cos(t)+sin(t);
k0=[0.0001 0.0001 0.0001 0.0001 0.0001 1.0009 0.00000001 10];
% x1 * sin(t) + x2 * (1 - 1 / (sqrt(1 - x3 * x3)) * exp(-x4) * sin(sqrt(1 -
% x3 * x3)*t + acos(x3)))++k(7)*cdf('poiss',t,k(8))
F=@(k,t)k(1) * sin(k(5)*t) + k(2) * (1 - 1 / (sqrt(1 - k(3) * k(3))) * exp(-k(4)) * sin(sqrt(1 - k(3) * k(3))*k(6)*t + acos(k(3))));
%% ==========lsqcurvefit==================

options = optimset('MaxFunEvals',20000);
[p,resnorm,residual,exitflag,output,lambda,jacobian]=lsqcurvefit(F,k0,t,Ca,[],[],options)
ci = nlparci(p,residual,jacobian)
yp=F(p,t);
fprintf('\n\n使用函数lsqcurvefit()估计得到的参数值为:\n')
fprintf('\t待求参数 k1 = %.6e\n',p(1))
fprintf('\t待求参数 k2 = %.6e\n',p(2))
fprintf('\t待求参数 k3 = %.6e\n',p(3))
fprintf('\t待求参数 k4 = %.6e\n',p(4))
fprintf('\t待求参数 k5 = %.6e\n',p(5))
fprintf('\t待求参数 k6 = %.6e\n',p(6))
% fprintf('\t待求参数 k7 = %.6e\n',p(7))
% fprintf('\t待求参数 k8 = %.6e\n',p(8))
fprintf('The sum of the squares is: %.3e\n\n',resnorm)
figure;plot(t,Ca,'r*-',t,yp(1,:),'bo-')
legend('exp Ca','pre Ca')

问题

参数越多，拟合效果就很不好