matlab练习程序（非线性常微分方程组参数拟合）

clear all;close all;clc;

bsrc = rand(4,1);
[t,h] = ode45(@(t,x)test1(t,x,bsrc),[0 100],[1 1]);
plot(t,h(:,1),‘r‘);
hold on;
plot(t,h(:,2),‘r‘);

hold on;
t = t(1:3:2*end/3);
x1 = h(1:3:2*end/3,1);
x2 = h(1:3:2*end/3,2);
plot(t,x1,‘ro‘);
plot(t,x2,‘ro‘);

T=min(t):0.1:max(t);        %插值处理，如果数据多，也可以不插值
X1=spline(t,x1,T)‘;
X2=spline(t,x2,T)‘;
plot(T,X1,‘b.‘);
plot(T,X2,‘b.‘);

dX1 = diff(X1)*10; dX1=[dX1;dX1(end)];
dX2 = diff(X2)*10; dX2=[dX2;dX2(end)];

%BFGS求解
syms k1 k2 k3 k4;
ff = sum((dX1 - (k1*exp(1./X1)+k4*sin(X2))).^2+(dX2-(cos(X1)*k2+k3*(1./X2))).^2);
dff1 = diff(ff,k1);dff2 = diff(ff,k2);
dff3 = diff(ff,k3);dff4 = diff(ff,k4);

f = inline(char(ff),‘k1‘,‘k2‘,‘k3‘,‘k4‘);
g1 = inline(char(dff1),‘k1‘,‘k2‘,‘k3‘,‘k4‘);
g2 = inline(char(dff2),‘k1‘,‘k2‘,‘k3‘,‘k4‘);
g3 = inline(char(dff3),‘k1‘,‘k2‘,‘k3‘,‘k4‘);
g4 = inline(char(dff4),‘k1‘,‘k2‘,‘k3‘,‘k4‘);

b = rand(4,1);
rho=0.2;sigma=0.2;
H=eye(4);

re=[];
for i=1:1000
    g=[g1(b(1),b(2),b(3),b(4));
        g2(b(1),b(2),b(3),b(4));
        g3(b(1),b(2),b(3),b(4));
        g4(b(1),b(2),b(3),b(4));];
    
    dk=-inv(H)*g;
    mk=1;
    
    for j=1:50
        new=f(b(1)+rho^j*dk(1),...
            b(2)+rho^j*dk(2),...
            b(3)+rho^j*dk(3),...
            b(4)+rho^j*dk(4));
        
        old=f(b(1),b(2),b(3),b(4));
        if new < old +sigma*rho^j*g‘*dk
            mk=j;
            break;
        end
    end
    
    norm(g)
    if norm(g)<1e-30 || isnan(new)
        break;
    end
    
    b = b + rho^mk*dk;
    gg=[g1(b(1),b(2),b(3),b(4));
        g2(b(1),b(2),b(3),b(4));
        g3(b(1),b(2),b(3),b(4));
        g4(b(1),b(2),b(3),b(4));];
    
    q = gg - g;             %q = g(k+1)-g(k)
    p = rho^mk*dk;          %p = x(k+1)-x(k)
    H = H - (H*p*p‘*H)/(p‘*H*p) + (q*q‘)/(q‘*p);    %矩阵更新
    
end

bsrc
b

[t,h] = ode45(@(t,x)test1(t,x,b),[0 100],[1 1]);
plot(t,h(:,1),‘g‘);
hold on;
plot(t,h(:,2),‘g‘);

function dy=test1(t,x,A)
a = A(1);       
b = A(2);     
c = A(3);
d = A(4);
dy=[a*exp(1/x(1))+d*sin(x(2));
   cos(x(1))*b+c/x(2)];