精华内容
下载资源
问答
  • 原问题地址:matlab 分式拟合,matlab 微分方程组的参数拟合_催眠神兔的博客-CSDN博客 微分方程式: x'=dx/dt=a*0.0321*(b-x)-d*x-dy/dt, y'= dy/dt=0.25*p1*exp(-p1*t)*x , 四个待求参数:a、b、d、p1 t、x、y...

    原问题地址:matlab 分式拟合,matlab 微分方程组的参数拟合_催眠神兔的博客-CSDN博客

    微分方程式:

    x'=dx/dt=a*0.0321*(b-x)-d*x-dy/dt,

    y'= dy/dt=0.25*p1*exp(-p1*t)*x ,

    四个待求参数:a、b、d、p1

    t、x、y数据见下面:

    //0 0 0  //这是初值

    0,0,0,

    0.1,0.486966799,0.048018378,

    0.167,1.6657,0.05823,

    0.2,0.860306078,0.060834243,

    0.3,1.156255213,0.064254733,

    0.4,1.390856542,0.065167644,

    0.5,1.67518,0.06638,

    0.6,1.724247244,0.065476325,

    0.7,1.841108525,0.065493681,

    0.8,1.93374543,0.065498314,

    0.9,2.007179471,0.06549955,

    1,1.92438,0.05641,

    1.1,2.111536196,0.065499968,

    1.2,2.148115682,0.065499991,

    1.3,2.177112544,0.065499998,

    1.4,2.200098596,0.065499999,

    1.5,2.218319829,0.0655,

    1.6,2.232763952,0.0655,

    1.7,2.244213928,0.0655,

    1.8,2.253290418,0.0655,

    1.9,2.260485427,0.0655,

    2,2.16156,0.07359,

    2.1,2.270710216,0.0655,

    2.2,2.274294246,0.0655,

    2.3,2.277135335,0.0655,

    2.4,2.27938749,0.0655,

    2.5,2.281172792,0.0655,

    2.6,2.282588016,0.0655,

    2.7,2.283709875,0.0655,

    2.8,2.284599183,0.0655,

    2.9,2.285304144,0.0655,

    3,2.44976,0.15449,

    3.1,2.286305961,0.0655,

    3.2,2.286657121,0.0655,

    3.3,2.286935489,0.0655,

    3.4,2.287156153,0.0655,

    3.5,2.287331076,0.0655,

    3.6,2.287469738,0.0655,

    3.7,2.287579657,0.0655,

    3.8,2.287666791,0.0655,

    3.9,2.287735862,0.0655,

    4,2.287790616,0.0655,

    4.1,2.287834019,0.0655,

    4.2,2.287868426,0.0655,

    4.3,2.2878957,0.0655,

    4.4,2.287917321,0.0655,

    4.5,2.287934459,0.0655,

    4.6,2.287948045,0.0655,

    4.7,2.287958815,0.0655,

    4.8,2.287967352,0.0655

    分析:

    微分方程中

    dx=a*0.0321*(b-x)-d*x-dy=0.0321*(a*b-a*x)-d*x-dy=0.0321*a*b-0.0321*a*x-d*x-dy=0.0321*ab-(0.0321*a-d)*x-dy=0.0321*ab-ad*x-dy

    等效于:
    dx=0.0321*ab-ad*x-dy

    故拟合参数a、b、d、p1中,前三个参数a、b、d不是独立的,如果进行拟合,将有无穷多最优解,故为过拟合现象。因而简化为三个拟合参数ab、ad、p1。

    Lu脚本代码(使用OpenLu+lugslmath计算,可从www.forcal.net下载):

    !!!using["luopt","math"]; //使用命名空间
    f(t,x,y,dx,dy, params :: ab,ad,p1)=
    {
        dy=0.25*p1*exp(-p1*t)*x,
        //dx=a*0.0321*(b-x)-d*x-dy, //dx=0.0321*(a*b-a*x)-d*x-dy=0.0321*a*b-0.0321*a*x-d*x-dy=0.0321*ab-(0.0321*a-d)*x-dy=0.0321*ab-ad*x-dy
        dx=0.0321*ab-ad*x-dy,
        0 //必须返回0
    };
    目标函数(_ab,_ad,_p1 : i,s,tyz : tyArray,tA,max, ab,ad,p1)=
    {
        ab=_ab, ad=_ad, p1=_p1, //传递优化变量
        //最后一个参数100表示gsl_ode函数在计算时,最多循环计算100次,这样可以提高速度
        tyz=gsl_ode[@f,nil,0.0,tA,ra1(0,0), 1e-6, 1e-6, gsl_rk8pd, 1e-6,100],
        i=0, s=0, while{++i<max, s=s+[tyz(i,1)-tyArray(i,1)]^2+[tyz(i,2)-tyArray(i,2)]^2},
        s
    };
    main(::tyArray,tA,max)=
    {
        tyArray=matrix{ //存放实验数据t、x、y
            "0,0,0,
    0.1,0.486966799,0.048018378,
    .........省略数据
    4.8,2.287967352,0.0655"
        },
        len[tyArray,0,&max], tA=tyArray(all:0), //用len函数取矩阵的行数,tA取矩阵的列
        Opt1[@目标函数] //Opt1函数全局优化
        //Opt[@目标函数] //也可以使用此语句
    };

    结果(ab,ad,p1,最小值):

    196.3535183967741         2.775743871883714         20.18977085052178         0.9084940395557356

    以下Lu代码可绘制图形:

    !!!using["luopt","math","win"]; //使用命名空间
    f(t,x,y,dx,dy, params :: ab,ad,p1)=
    {
        dy=0.25*p1*exp(-p1*t)*x,
        //dx=a*0.0321*(b-x)-d*x-dy, //dx=0.0321*(a*b-a*x)-d*x-dy=0.0321*a*b-0.0321*a*x-d*x-dy=0.0321*ab-(0.0321*a-d)*x-dy=0.0321*ab-ad*x-dy
        dx=0.0321*ab-ad*x-dy,
        0 //必须返回0
    };
    目标函数(_ab,_ad,_p1 : i,s,tyz : tyArray,tA,max, ab,ad,p1)=
    {
        ab=_ab, ad=_ad, p1=_p1, //传递优化变量
        //最后一个参数100表示gsl_ode函数在计算时,最多循环计算100次,这样可以提高速度
        tyz=gsl_ode[@f,nil,0.0,tA,ra1(0,0), 1e-6, 1e-6, gsl_rk8pd, 1e-6,100],
        i=0, s=0, while{++i<max, s=s+[tyz(i,1)-tyArray(i,1)]^2+[tyz(i,2)-tyArray(i,2)]^2},
        s
    };
    init(main:tyz:tyArray,tA,max,ab,ad,p1)=
    {
        tyArray=matrix{ //存放实验数据t、x、y
            "0,0,0,
    0.1,0.486966799,0.048018378,
    
    ... ...省略数据
    
    4.8,2.287967352,0.0655"
        },
    
        len[tyArray,0,&max], tA=tyArray(all:0), //用len函数取矩阵的行数,tA取矩阵的列
        ab=0.0, ad=0.0, p1=0.0, Opt1[@目标函数, optstarter,&ab,&ad,&p1,0],   //Opt1函数全局优化
        tyz=gsl_ode[@f,nil,0.0,tA,ra1(0,0), 1e-6, 1e-6, 2, 1e-6,100],  //获得最终结果
        outa[tyz],                  //outa输出结果
        luShareX2[tyArray, tyz]     //绘制共享X轴视图
    };
    ChartWnd[@init];                //显示窗口并初始化

    图形:

    展开全文
  • function k1k2k32 format long clear all clc tspan = [0 6 24 44 68 72 74 92 104 116]';%%这是时间 yexp= [3.111,3.639,3....以及想用fitnlm函数进行拟合,达到一些参数评价t值p值等,但用fitnlm运行不起来。 @beefly

    function k1k2k32

    format long

    clear all

    clc

    tspan = [0 6 24 44 68 72 74 92 104 116]';%%这是时间

    yexp= [3.111,3.639,3.887,4.289,4.658,5.531,6.218,6.979,7.111,7.114]';%%%这是菌落总数

    x0 = [3.111];

    k0 = [0.1  3  45 4 37 8.1];

    lb = [0  0  0  0 0  0];

    ub = [100  100  100 100 100 100];

    %

    % % opts = statset('nlinfit');

    % % opts.RobustWgtFun = 'bisquare';

    % % mdll = fitnlm(tspan,yexp,ObjFunc,b0,'Options')

    % mdll = fitnlm(tspan,yexp,ObjFunc,b0,'Options',opts)

    options = optimoptions(@lsqnonlin,'Algorithm','trust-region-reflective');

    [k,resnorm,residual,exitflag,output,lambda,jacobian] = ...

    lsqnonlin(@ObjFunc,k0,lb,ub,[],tspan,x0,yexp);

    ci = nlparci(k,residual,jacobian);

    fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')

    fprintf('\tk1 = %.9f ± %.9f\n',k(1),ci(1,2)-k(1))

    fprintf('\tk2 = %.9f ± %.9f\n',k(2),ci(2,2)-k(2))

    fprintf('\tk3 = %.9f ± %.9f\n',k(3),ci(3,2)-k(3))

    fprintf('\tk4 = %.9f ± %.9f\n',k(4),ci(4,2)-k(4))

    fprintf('\tk5 = %.9f ± %.9f\n',k(5),ci(5,2)-k(5))

    fprintf('\tk6 = %.9f ± %.9f\n',k(6),ci(6,2)-k(6))

    fprintf('  The sum of the squares is: %.9e\n\n',resnorm)

    tsa=0:0.01:max(tspan);

    [tsa ysa]=ode45(@KineticsEqs,tsa,x0,[],k);

    figure(1),

    plot(tsa,ysa(:,1),'b',tspan,yexp,'or'),legend('计算值','实验值','Location','best');

    function f = ObjFunc(k,tspan,x0,yexp)           % 目标函数

    [tspan  y] = ode45(@KineticsEqs,tspan,x0,[],k);

    f = (y-yexp);

    function dCdt = KineticsEqs(tspan ,C,k)

    T=4*(tspan >=0 & tspan <=30)+10*(tspan >31 & tspan <=90)+4*(tspan >90 ); % 假设这是我写的温度关于时间的分段函数

    umax=((k(2)*(T-k(3))*(T-k(4))*(T-k(4)))/(((k(5)-k(4))*(T-k(5))-(k(5)-k(3))*(k(5)+k(4)-2*T))*(k(5)-k(4))));

    dCdt=(1/(1+exp(-4*(tspan -k(1)))))*umax*(1-exp(C-k(6)));

    这是根据他人的代码更改的,得到的结果是置信区间过大?想有什么办法解决?以及想用fitnlm函数进行拟合,达到一些参数评价t值p值等,但用fitnlm运行不起来。 @beefly

    展开全文
  • 如何使用Matlab编程进行参数拟合 1前言2基本概念和原理3主要内容4实例5涉及的文件1前言之前帮疯学网做过一个利用Matlab编程进行参数拟合 的教程,由于疯学网好像倒闭了,希望之前做的工作不要白费,这里...

    如何使用Matlab编程进行参数拟合 

     

    1前言
    2基本概念和原理
    3主要内容
    4实例
    5涉及的文件



    1前言


    之前帮疯学网做过一个利用Matlab编程进行参数拟合 的教程,由于疯学网好像倒闭了,希望之前做的工作不要白费,这里拿出来分享下,希望能对虫友的学习、科研工作有所帮助。其他的不多说,言归正传,下面从原理和实例对如何使用Matlab编程进行参数拟合进行讲解。


    2基本概念和原理

    所谓参数拟合,就是已知试验或者真实数据,然后寻找一个模型对其规律进行模拟,求取模型中未知参数的一个过程。根据模型的复杂程度我分成以下几类讲解:
    代数方程模型
    线性模型
    非线性模型
    单变量,单目标问题
    多变量,单目标问题
    多变量,多目标问题,共享参数
    复数模型拟合
    微分方程模型
    简单微分方程参数拟合
    复杂微分方程参数拟合

    3主要内容



    3.1代数方程模型

    y=f(x,β)
    x, n维自变量x=[x1,x2,…,xn]'
    β,p维参数向量β =[β1, β2,…, βn]'
    y,m维因变量y=[y1,y2,…,yn]'
    f,m维函数向量f=[f1,f2,…,fn]'

    Matlab 求解函数
    线性最小二乘法:ployfit, regress
    非线性最小二乘法:fit, lsqnonlin,lsqcurvefit,nlinfit
    nlintool,fmincon


    3.2微分方程模型

    dx/dt=f(t,x,β,u)  x(t0)=x0
    x, n维状态变量x=[x1,x2,…,xn]'
    x0, n维状态变量初值x=[x10,x20,…,xn0]'
    β,p维参数向量β =[β1, β2,…, βp]'
    u,r维操作变量y=[u1,u2,…,um]'
    f,ODE方程m维函数向量f=[f1,f2,…,fm]'

    对于给定β,由龙格-库塔积分可得x(t)
    y(t)=g(x(t), β) 
    y为m维输出量y=[y1,y2,…,yn]'
    g为输出量y与状态变量x之间非线性关系的函数向量g=[g1,g2,…,gn]'
    用导数法和积分法将ODE问题转化为代数方程问题


    3.2优化准则和参数初值确定方法

    优化准则:最小二乘法,极大似然,概率密度函数
    非线性模型必须采用非线性回归的方法
    参数初值确定方法
    (1)根据模型结构 和本质
    描述物理系统的模型的结构和本质可能指明未知参数的取值范围
    (2)模型方程的渐进行为
    例如指数衰减模型yi=k1+k2exp(-k3*xi)
    xi-->∞,yi约为k1, xi-->0,yi=k1+k2(1-k3*xi),利用接近0的x及对应y回归,结合k1,就可得到k1 k2 k3初值
    (3)模型方程的变换
    把非线性系统转化为线性系统,线性最小二乘结果作为非线性估计的初值
    (4)直接搜索,全局算法
    如果对参数不了解,可用直接搜索或者全局(多起点,遗传等算法处理)


    3.4决定性指标



    回归均值的总偏离 
    Ne:实验点数目,
    yei,yci,i次实验值与计算值

    由于公式比较难显示,见附件ppt里面内容
    如何使用Matlab编程进行参数拟合

    4实例



    4.1线性拟合函数:regress()

    调用格式:b=regress(y,X)
    [b,bint,r,rint,stats]= regress(y,X)
    [b,bint,r,rint,stats]= regress(y,X,alpha)
    该函数求解线性模型:y=Xβ+ε
    β是p1的参数向量;ε是服从标准正态分布的随机干扰的n1的向量;y为n1的因变量向量;X为np自变量矩阵。
    bint返回β的95%的置信区间。r中为形状残差,rint中返回每一个残差的95%置信区间。Stats向量包含R2统计量、回归的F值和p值。

    例 设y的值为给定的x的线性函数加服从标准正态分布的随机干扰值得到。即y=10+x+ε ;求线性拟合方程系数。
    x=[ones(10,1) (1:10)']
    y=x*[10;1]+normrnd(0,0.1,10,1)
    [b,bint, r,rint,stats]=regress(y,x,0.05)
    rcoplot(r,rint)


    4.2简单线性模型-多项式拟合


    多项式曲线拟合函数:polyfit( )
    调用格式:p=polyfit(x,y,n)
                    [p,s]= polyfit(x,y,n)
    x,y为数据点,n为多项式阶数,返回p为幂次从高到低的多项式系数向量p。矩阵s用于生成预测值的误差估计。
    多项式曲线求值函数:polyval( )
    调用格式:y=polyval(p,x)
                    [y,DELTA]=polyval(p,x,s)
    y=polyval(p,x)为返回对应自变量x在给定系数P的多项式的值
    多项式曲线拟合的评价和置信区间函数:polyconf( )
    调用格式:[Y,DELTA]=polyconf(p,x,s)
                    [Y,DELTA]=polyconf(p,x,s,alpha)
    多项式输出 ps = poly2str(p,'x')    
    codefile Fit_polynomial  
                          

    4.3稳健回归函数:robustfit( )


    稳健回归是指此回归方法相对于其他回归方法而言,受异常值的影响较小。
    调用格式:        
    b=robustfit(x,y)
    [b,stats]=robustfit(x,y)
    [b,stats]=robustfit(x,y,’wfun’,tune,’const’)

    例:演示一个异常数据点如何影响最小二乘拟合值与稳健拟合。首先利用函数y=10-2x加上一些随机干扰的项生成数据集,然后改变一个y的值形成异常值。调用不同的拟合函数,通过图形观查影响程度。
    code File Robust_Fit


    4.4 非线性问题使用matlab函数-lsqcurvefit,lsqnonlin


    [x resnorm] = lsqcurvefit(fun,x0,xdata,ydata,...)
    fun   是我们需要拟合的函数,
    x0    是我们对函数中各参数的猜想值,
    xdata  则是自变量的值
    ydata  是因变量的值,目标值
    min  sum {(FUN(X,XDATA)-YDATA).^2} 

    x = lsqnonlin(fun,x0)   
     x0为初始解向量;fun为优化函数,fun返回向量值F,而不是平方和值,平方和隐含在算法中,fun的定义与前面相。 
    x = lsqnonlin(fun,x0,lb,ub)     
    lb、ub定义x的下界和上界
    x = lsqnonlin(fun,x0,lb,ub,options)    
    options为指定优化参数,若x没有界,则lb=[ ],ub=[ ]
    min  sum {FUN(X).^2}


    4.5非线性问题使用matlab函数-fmincon


    x= fmincon(fun,x0,A,b) 
    给定初值x0,求解fun函数的最小值x。fun函数的约束条件为A*x<= b,x0可以是标量或向量。x= fmincon(fun,x0,A,b,Aeq,beq) 
    最小化fun函数,约束条件为Aeq*x= beq 和 A*x <= b。若没有不等式线性约束存在,则设置A=[]、b=[]。x= fmincon(fun,x0,A,b,Aeq,beq,lb,ub) 
    定义设计变量x的线性不等式约束下界lb和上界ub,使得总是有lb<= x <= ub。若无等式线性约束存在,则令Aeq=[]、beq=[]。x= fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon)
    在上面的基础上,在nonlcon参数中提供非线性不等式c(x)或等式ceq(x)。fmincon函数要求c(x) <= 0且ceq(x)= 0。x= fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options) 
    用options参数指定的参数进行最小化。x= fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options,P1,P2,...) 将问题参数P1, P2等直接传递给函数fun和nonlin。若不需要这些变量,则传递空矩阵到A, b, Aeq, beq, lb, ub, nonlcon和 options
    min F(X) fun函数返回一向量或者标量


    4.6非线性问题线性化

    如何使用Matlab编程进行参数拟合-1


    4.7单变量,单目标问题(cftool 的应用)

    例子
    渗流公式为y=A*(x-xc)^p其中x为自变量,y为因变量,A、xc和p均为常数
    为了测试模拟,设定A=18.5,xc=0.095,p=-2.3,得到以下数据
    x                y-----------------------
    -0.1001        3.5E+06
    0.1002        3.3E+06
    0.11           2.9E+05
    0.12           9.0E+04
    0.15           1.5E+04
    0.2             3.3E+03
    0.3             7.1E+02
    0.4             2.8E+02
    0.5             1.5E+02
    0.6             8.9E+01
    http://muchong.com/bbs/viewthread.php?tid=3866180
    code file  percolation_fit


    4.8 复数模型拟合

    将模型分离成实部和虚部
    求解如下优化模型

    例子
    http://muchong.com/bbs/viewthread.php?tid=4268659&target=self&page=1
    模型
    变量参数为:a,b,c,d1,m1,n1,d2,m2,n2,d3,m3,n3,拟合曲线为 e=e1+ie2=a-b/(x^2+icx)+m1/(n1-x^2-id1x)+m2/(n2-x^2-id2x)+m3/(n3-x^2-id3x)
    Data
        x                e1             e2
    5.16957        1.854        4.5139
    4.96279        1.9351        4.5777
    4.77192        2.0221        4.4781

    code file  complexfit


    4.9简单微分方程参数拟合


    1.由给定的ODE参数初值结合ODE方程dC/dt=f(k,c,t),利用ode45积分可得对应时间点的浓度预测值Cp
    2.以浓度的预测值Cp与实验值Ce的残差平方和为优化目标,利用lsqnonlin或者fmincon进行搜索获得最优的ODE参数,每搜一次,调用ode45计算预测浓度值,直到优化完成


    http://muchong.com/bbs/viewthread.php?tid=4685934&target=self&page=1
    问题
    实验数据为:(t,c)=(0,0.69)(2,0.645)(4,0.635)(8,0.62)(24,0.61)(48,0.61).
    其中t为时间,c为某离子的浓度。动力学方程模型为:-dc/dt=k*(c0-c)^(1/3)*(c-c~).其中c0为初始浓度可以取0.7,c~为平衡浓度取0.61.
    code file ODE_parafit


    4.10复杂微分方程参数拟合


    在一个等温间歇反应器中进行复杂反应,描述该反应系统的七个微分方程
    有5个参数k1, k2, k3, k4, k5, 初始状态x(0)=[0.1883 0.2507 0.0467 0.0899 0.1804 0.1394 0.1046]',根据下表数据用最优化方法估计参数k
           t            y1(x1)    y2(x4)    y3(x5)    y4(x6)  
           0        0.1883    0.0899    0.1804    0.1394
        0.0100    0.2047    0.0866    0.1729    0.1297
        0.0200    0.2181    0.0856    0.1680    0.1205
    …….
    微分方程:









    code file: KineticsEst5.m


    4.11 多元非线性拟合,全局优化


    例子
    https://zhidao.baidu.com/question/168077392.html
    matlab nlinfit多元非线性拟合 错在哪里?
    需要拟合如下形式的模型:
    y = beta(1)+beta(2)*V+exp(beta(3)+beta(4)*R)

    我使用的代码如下:
    >> V=[47,69,54,65.6,46,61,66.5,50.8,55.9,44.5,53.3,54,56.8,51.9,60.1,67.3,62.4,61.1,61.7,76,51.7,50.7,57.1,65.6,53.1,60.6,59.1,46,49,47.3,75.3,52.5,58.3,54.2,46.8,55.1,47.8];
    >> R=[47,1081,186,127,50,109,397,178,189,72,86,98,273,186,202,137,138,104,672,270,94,57,59,187,129,55,288,77,155,65,972,85,884,195,105,172,97];
    >> x = [V;R]';
    >> MSR=[18,11,11,16,20,23,26,2.9,11.1,12.2,22.1,19.8,8.6,4.3,17.6,35.9,27.2,30.7,11.9,42.3,16.7,31.9,41.6,28.2,12.2,50.9,12.2,12.6,1.9,20.7,34.6,21,5.1,7.7,5.4,10.9,9.2];
    >> y = MSR';
    beta=nlinfit(x,y,myfun,[-100,1,-10,1000]

    m-函数为:
    function y  = myfun(beta,x)
    y = beta(1)+beta(2)*x(:,1)+exp(beta(3)+beta(4)*x(:,2));

    错误的原因是:Input argument "beta" is undefined.beta没有定义?
    请问错在哪里?在线等待哦,谢谢!
    在命令窗口输入yy=myfun(beta,x),回车运行试试,
    也是不行的,

    ??? Error using ==> beta at 21
    Not enough input arguments.


    code file:Muti_var_fitCODE:
    function Muti_var_fit
    % matlab nlinfit多元非线性拟合错在哪里?
    % 举报|2010-07-19 11:59transtorseu | 分类:其他编程语言 | 浏览1500次
    % 需要拟合如下形式的模型:
    %y = beta(1)+beta(2)*V+exp(beta(3)+beta(4)*R)
    %http://zhidao.baidu.com/link?url=7Jue1nv1dhSE2OpGMv6pKWOW7NX0zgmrpAZV6usGpavPA8PSUJSWY0Hp0AdgTkdbmdBTnYnLaIJXg4Z8wt2r8_
    myfun=@(x,xdata)x(1)+x(2)*xdata(:,1)+exp(x(3)+x(4)*xdata(:,2));
    V=[47,69,54,65.6,46,61,66.5,50.8,55.9,44.5,53.3,54,56.8,51.9,60.1,67.3,62.4,61.1,61.7,76,51.7,50.7,57.1,65.6,53.1,60.6,59.1,46,49,47.3,75.3,52.5,58.3,54.2,46.8,55.1,47.8]';
    R=[47,1081,186,127,50,109,397,178,189,72,86,98,273,186,202,137,138,104,672,270,94,57,59,187,129,55,288,77,155,65,972,85,884,195,105,172,97]';
    xdata = [V R];
    MSR=[18,11,11,16,20,23,26,2.9,11.1,12.2,22.1,19.8,8.6,4.3,17.6,35.9,27.2,30.7,11.9,42.3,16.7,31.9,41.6,28.2,12.2,50.9,12.2,12.6,1.9,20.7,34.6,21,5.1,7.7,5.4,10.9,9.2];
    ydata = MSR';
    %beta0=[-80 1 4 -0.01]
    x0=[-20 1 1 0.01];
    % [x,residual,jacobian,covb,resnorm]=nlinfit(xdata,ydata,myfun,x0);
    % ci = nlparci(x,residual,'covar',covb)
    %% =======利用lsqcurvefit 非线性最小二乘法
    [x,resnorm,residual,exitflag,output,lambda,jacobian]=lsqcurvefit(myfun,x0,xdata,ydata);
    ci = nlparci(x,residual,jacobian)
    fprintf('\n\n使用函数lsqcurvefit()估计得到的参数值为:\n')
    fprintf('\t待求参数 k1 = %.6f\n',x(1))
    fprintf('\t待求参数 k2 = %.6f\n',x(2))
    fprintf('\t待求参数 k3 = %.6f\n',x(3))
    fprintf('\t待求参数 k4 = %.6f\n',x(4))
    fprintf('  The sum of the squares is: %.3e\n\n',resnorm)
    figure;plot(1:numel(ydata),ydata,'r-*',1:numel(ydata),myfun(x,xdata),'bo-')
    legend('real','model')
    Text=['使用局部优化函数lsqcurvefit估计得到的参数'];
    title(Text)
    %% ==========利用globalsearch和 fmincon=========
    tic
    x0=[-10 1 1 0.1];
    F =@(x)norm(x(1)+x(2)*xdata(:,1)+exp(x(3)+x(4)*xdata(:,2))-ydata);
    gs = GlobalSearch('Display','iter');
    opts=optimset('Algorithm','trust-region-reflective');
    problem=createOptimProblem('fmincon','x0',x0,...
        'objective',F,'lb',[-100,-100,-100,-1],'ub',[100,100,100,1],'options',opts);
    [xming,fming,flagg,outptg,manyminsg] = run(gs,problem)
    t1=toc
    figure;plot(1:numel(ydata),ydata,'r-*',1:numel(ydata),myfun(xming,xdata),'bo-')
    legend('real','model')
    Text1=['利用全局算法globalsearch和 fmincon估计得到的参数,耗时',num2str(t1),'s'];
    title(Text1)
    %% =========全局算法 multistart和lsqcurvefit
    tic
    ms=MultiStart;
    opts=optimset('Algorithm', 'trust-region-reflective');
    problem=createOptimProblem('lsqcurvefit','x0',x0,'xdata',...
        xdata, 'ydata',ydata,'objective',myfun,'lb',[-100,-100,-100,-1],'ub',[100,100,100,1],'options',opts);
    [xminm,fminm,flagm,outptm,manyminsm]=run(ms,problem,20)
    t2=toc
    figure;plot(1:numel(ydata),ydata,'r-*',1:numel(ydata),myfun(xminm,xdata),'bo-')
    legend('real','model')
    Text2=['利用全局算法 multistart和lsqcurvefit估计得到的参数,耗时',num2str(t2),'s'];
    title(Text2)
    %% =============遗传算法的参数识别=======
    tic
    Fgen =@(x)norm(x(1)+x(2)*xdata(:,1)+exp(x(3)+x(4)*xdata(:,2))-ydata);
    options = gaoptimset('TolFun',1e-10);
    [xgen,fval] = ga(Fgen,4,[],[],[],[],[-100;-100;-100;-1], [100;100;100;1])
    [xgen,resnorm,residual,exitflag,output,lambda,jacobian]=lsqcurvefit(myfun,xgen,xdata,ydata);
    ci = nlparci(xgen,residual,jacobian)
    t3=toc
    figure;plot(1:numel(ydata),ydata,'r-*',1:numel(ydata),myfun(xgen,xdata),'bo-')
    legend('real','model')
    Text3=['利用全局算法遗传算法的参数识别估计得到的参数,耗时',num2str(t3),'s'];
    title(Text3)

    结果如下
    如何使用Matlab编程进行参数拟合-2
    如何使用Matlab编程进行参数拟合-3
    如何使用Matlab编程进行参数拟合-4
    如何使用Matlab编程进行参数拟合-5
    由图可全局算法估计参数结果优于一次非线性优化估计的参数


    5涉及的文件 

    (见附件)本贴录有视频,如需要私信。

    转载于:https://www.cnblogs.com/Simulation-Campus/p/9001499.html

    展开全文
  • 第五章 控制系统仿真 5.2 微分方程求解方法 以一个自由振动系统实例为例进行讨论如下图 1 所示弹簧 - 阻尼系统参数如下 M=5 kg, b=1 N.s/m, k=2 N/m, F=1N x b M F k 图 1 弹簧-阻尼系统 假设初始条件为 t0 0 时将 m...
  • 参数的常微分方程拟合问题研究一.问题的背景:二.提出一个较为简单,但是很有代表性的一个问题:三.求解的基本原理:四.求解的基本算法:1.利用matlab遗传算法求解出其最优误差函数值:2.差分后用LINGO求解: 一....

    一.问题的背景:

    一般我们平常学习到的拟合问题大部分都是已知了函数表达式后,对于函数表达式里的系数进行的拟合问题。大部分的教材上其实并没有见讲解该如何对不能求解出相关表达式的函数关系的参数拟合求解。而且LINGO也是没有现成的微分符号表达2阶微分。因此,我通过在网上查阅部分资料提出了自己的一些小想法,用来求解一般的较为复杂的微分方程拟合的问题,并且使得解有一定的精度。

    二.提出一个较为简单,但是很有代表性的一个问题:

    d 2 y d t 2 = a y , 其 中 y ( 1 ) = 0.909 , y ( 1.1 ) = 0.706 \frac{d^2y}{dt^2}=ay, 其中y(1) = 0.909,y(1.1) = 0.706 dt2d2y=ay,y(1)=0.909,y(1.1)=0.706

    ​ 其中, t t t每隔 0.1 0.1 0.1个长度的数据如下:
    y r e a l = [ 0909 , 0.706 , 0.363 , 0.079 , − 0.226 , − 0.474 , − 0.765 , − 1.183 , − 1.436 , − 1.573 ] y_{real} = [0909,0.706,0.363,0.079,-0.226,-0.474,-0.765,-1.183,-1.436,-1.573] yreal=[0909,0.706,0.363,0.079,0.226,0.474,0.765,1.183,1.436,1.573]
    ​ 求: a ^ \hat{a} a^

    ​ 解决这类问题首先是不能用 m a t l a b matlab matlab o d e ode ode函数求解的,因为这是一个有 2 2 2点初值的问题,在 t = 1 t=1 t=1时刻 y y y的导数是不知道的。

    三.求解的基本原理:

    ​ 求解这个问题,而且包括更广放的一般的拟合问题,都要用到我们伟大的高斯老师的最小二乘法作为基本的工具。
    m i n a ∈ R ∑ i = 1 n ( y c a l c ( i ) − y r e a l ( i ) ) 2 min_{a\in R}\sum_{i=1}^n(y_{calc}^{(i)}-y_{real}^{(i)})^2 minaRi=1n(ycalc(i)yreal(i))2
    ​ 然而,要想求出计算值来实际上得先把这个微分问题转化为差分方程问题如下图:
    h = 0.1 ; t 0 = 1 ; y t 0 = 0.909 , y t 0 + h = 0.706 y t + h + y t − h − 2 y t h 2 = a y t h = 0.1;t_0 = 1;y_{t_0} = 0.909,y_{t_0+h} = 0.706\\ \frac{y_{t+h}+y_{t-h}-2y_t}{h^2} = ay_t h=0.1;t0=1;yt0=0.909,yt0+h=0.706h2yt+h+yth2yt=ayt
    ​ 这样只要我们知道了 a a a的大小就可以把每一点的 y y y值算出来了。这样其实就转化成了一个求 a a a值的无约束单目标优化问题, 目标函数就是我们的最小二乘函数

    四.求解的基本算法:

    1.利用 m a t l a b matlab matlab遗传算法求解出其最优误差函数值:

    ​ 我们利用 m a t l a b matlab matlab自带的现成的 g a ga ga工具箱求解,基本的代码如下:

    function err = obj_ode(a)
    h = 0.1;
    fun_ode = zeros(1,10);
    real_value = [0.909 0.706 0.363 0.079 -0.226 -0.474 -0.765 -1.183 -1.436 -1.573];
    fun_ode(1) = 0.909;
    fun_ode(2) = 0.706;
    for t = 3:10
        fun_ode(t) = (a*h^2+2)*fun_ode(t-1)-fun_ode(t-2);
    end
    err = sum((fun_ode-real_value).^2);
    end
    

    ​ 上面的函数是为了求出 a a a值固定时候的误差大小,也就是我们的最小二乘函数,现在要使得这个函数最小,就应该在主函数里用遗传算法来求解:

    clc,clear;
    real_value = [0.909 0.706 0.363 0.079 -0.226 -0.474 -0.765 -1.183 -1.436 -1.573];
    [a,error_value] = ga(@obj_ode,1);
    f = @(t,u)[u(2);a*u(1)^2];
    y_init = [0.909;(0.706-0.909)/0.1];
    [u,y] = ode45(f,1:0.1:1.9,y_init);
    plot(1:0.1:1.9,y(:,1),'r-',1:0.1:1.9,real_value,'b-');
    xlabel('自变量t');
    ylabel('因变量y');
    legend('计算值','真实值');
    fprintf('a的值是:%f',a);
    fprintf('最小误差的值是:%f',error_value);
    

    ​ 求出来最后的拟合曲线如下:

    在这里插入图片描述

    ​ 求出来的 a a a − 8.13 -8.13 8.13,误差在 0.49 0.49 0.49左右。但是从曲线上看貌似到最后误差变得有点大了

    2.差分后用 L I N G O LINGO LINGO求解:

    model:
    sets:
    time/1..100/:t,y;!100等分点;
    num/1..10/:y_real,y_calc;
    endsets
    data:
    y_real=0.909 0.706 0.363 0.079 -0.226 -0.474 -0.765 -1.183 -1.436 -1.573;!实际值;
    enddata
    CALC:
    h = 0.01;!步长;
    y(1) = 0.909;
    y(11) = 0.706;
    endcalc
    min = @sum(num:(y_real-y_calc)^2);
    @for(time(i)|i#gt#2:y(i)=(a*h^2+2)*y(i-1)-y(i-2));!差分关系;
    @for(num(i):y_calc(i) = y(10*(i-1)+1));!映射关系;
    @free(a);
    end
    

    ​ 最后求出的结果如下:

    在这里插入图片描述

    ​ 求出的 a a a 8.43 8.43 8.43,很明显示 L I N G O LINGO LINGO收敛到了局部最优解。没有得到我们想要的全局最优解。在这方面,遗传算法要比 L I N G O LINGO LINGO强上一些。

    其实也可以尝试以下用4阶龙格库塔法编写 L I N G O LINGO LINGO求解,有尝试过的小伙伴可以自己比较一下结果哦!

    展开全文
  • R^2): 0.572359944199415 确定系数(DC): 0.536839307201338 F统计(F-Statistic): 160.74113587027 参数 最佳估算 -------------------- ------------- a 11.4149902398646 b -11.2943273492988 c 38....
  • 根据两个式子(一个隐函数式微分方程,一个代数方程)拟合7个未知参数fac1~fac7。式中lambv'=d(lambv)/d(lamb),lambv为中间变量。lambv'=((1.0/3/fac7)*(fac1*((lamb/lambv')^fac2-(lambv'/lamb)^(0.5*fac2))+fac3*...
  • 这是我所要用的二阶微分方程,现已通过试验得到数据,想通过该方程拟合数据得到方程中的四个参数,通过查询有说二阶微分方程转化为差分方程或者积分方程,然后MATLAB编程拟合,想知道如何转化为差分方程或者积分方程...
  • x,y,z关于t的微分方程组中存在未知参数,已知x,y,z,t的多组离散数据,能求出方程组中各个参数值么
  • MATLAB使用欧拉Euler法求解微分方程组 部分源码 clear;clc c=2/3; %设置c的值 x(1)=0.1; %设置x初值为0.1 y(1)=0.3; %设置y初值为0.3 h=0.05; %设置步长为0.05
  • 参数微分方程matlab求解

    千次阅读 2018-04-13 13:53:54
    参数微分方程matlab求解 对于Lorenz微分方程组 ⎧⎩⎨⎪⎪x˙1=−83x1+x2x3x˙2=−10x2+10x3x˙3=−x1x2+28x2−x3{x˙1=−83x1+x2x3x˙2=−10x2+10x3x˙3=−x1x2+28x2−x3\begin{cases} \dot{x}_{1} = - \frac{8...
  • matlab微分方程代码ESO208编程分配 此仓库维护2019-20 I学期ESO208A(计算方法和数值分析)作业的代码。ESO208A总共包括4个编程作业。 这些程序已在MATLAB版本R2019a上编写。 每个PA文件夹均包含分配问题说明,...
  • 中山大学 数学实验与数学软件 第08章 MATLAB微分方程数值解法(共33页).pptx 中山大学 数学实验与数学软件 第09章 MATLAB数值线性代数(共29页).pptx 中山大学 数学实验与数学软件 第10章 MATLAB进阶程序设计与...
  • 事实证明,显式微分方程求解器在拟合React工程动力学模型方面非常有用。 我确实相信持续改进,因此我同意我可以使这些 mfiles 更好。 请随时提供有用的意见。 这些 mfile 最多可以处理大约 10 个独立参数。 尽管...
  • 在数学建模大赛中,常微分方程的求解至关重要。以下是作者在学习过程中,以SARS病毒传染病SIR模型为例进行的数据处理以及模型求解的过程。
  • %改进欧拉法求解常微分方程数值解 %初始条件 %dy/dx=f(x) %y(x0)=y0 clc;clear; f=@(x,y)1/x+1/y;%方程 x=[];y=[];%初值 x(1)=1;y(1)=1; num=100;%迭代次数 a=1;b=50;%自变量区间 h=(b-a)/(num-1); for i=2:num x(i...
  • matlab求解微分方程组代码 计算方法(calculation methodb) 该项目是《计算方法》一书中提到的经典方法和算法的matlab程序实现,包含代码详解和运行过程。 :grinning_face_with_big_eyes: 1.简介 2.线性方程组的...
  • matlab求解常微分方程组/传染病模型并绘制SIR曲线

    万次阅读 多人点赞 2016-04-19 14:34:42
    本文用matlab求解常微分方程组。 看了很多关于传染病模型的matlab程序,大都是绘制出两条曲线(I、S)的,本文最大的不同是绘出SIR三条曲线。 先给出SIR微分方程组   函数文件: run的程序: matlab输出: ...
  • 有位小伙伴在matlab编程爱好者(群号:531421022)群中问道有关时滞微分方程matlab解法,问题是选自由清华大学出版社出版、薛定宇著的《高等应用数学问题的MATLAB求解 (第四版)》的课后习题,问题的如下:显然这是...
  • **前言:**大学期间只学习过《常微分方程》,没想到有些学校竟然还学《时滞微分方程》,于是找到一本由内藤敏机(日本)等著,...没找到有人用欧拉法解一阶时滞微分方程的,于是一不做二不休便用MATLAB实现了一下下...
  • 微分方程的形式为: dy/dx=f(x,y) y(x0)=y0 基本思想: 向前欧拉差分方法: (yn+1-yn)=hf(xn,yn) 向后欧拉差分方法: (yn-yn-1)=hf(xn,yn) 泰勒级数方法:
  • 问题来源:如何拟合微分方程组的参数? – MATLAB中文论坛 (ilovematlab.cn) 微分方程组如下: dx/dt=a*x-b*x*y dy/dt=-c*y+d*x*y 数据如下: t x(t) y(t) 11 45.79 41.40 12 53.03 38.90 13 64.05 36.78 14 ...
  • MATLAB与最小二乘法拟合数据

    千次阅读 2021-04-04 19:57:41
    1、前言:学习并记录的原因 ...(1)常见的数据拟合有直线拟合、多项式拟合、插值拟合等首先利用MATLAB中的函数来直观体验下一次拟合。 参考:http://www.qinms.com/work/nihe.html(拟合方法) h...
  • Matlab计算微分方程曲线求导及过曲线上点的切线方程求解f(x)=x^2一元二次方程上某点的切线方程并绘制出方程的切线图。点(4,f(4))是曲线方程f(x)上的一个点,求出该点的切线并绘制出来。画出f(x)= x^2方程曲线。对f(x...
  • MATLAB程序实现差分进化算法,程序包含5个文件,分别为主程序、初始化种群、适应度函数(选择)、交叉、变异。程序示例为设计一阶控制器使离散传递函数(z-1)(z+0.3)/z(z-2)(z-0.5)稳定。
  • 在数据处理中经常会需要对数据进行拟合拟合完成之后可以通过拟合曲线的方程对数据进行预测。下面主要介绍一下如何适用matlab自带的拟合工具包对数据进行拟合,全程不需要编写一句代码,拟合完成之后还能生成函数...
  • 在python里如何用4阶龙格库塔法解微分方程组呢?(在matlab里是写一个函数文件,然后ode45调用并赋初值) 大佬可否根据实例在python写一段程序我参考学习一下 (最后还需要绘制曲线t...
  • 微分方程建模详解

    千次阅读 2021-01-25 00:00:00
    微分方程建模
  • 更多内容尽在个人专栏:matlab学习之前我们学习过matlab的各种解方程的函数工具,这一节我们再来学习一种,(常)微分方程的求解工具。微积分的符号表示方法:我们先了解一下怎么用符号表示导数这个变量:符号计算中...

空空如也

空空如也

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

matlab微分方程参数拟合

matlab 订阅