matlab練習程序（非線性常微分方程組參數擬合）

和線性常微分方程組參數擬合類似，我們要用差分代替微分，然后進行插值處理，然后構造最小化函數。

最后用最優化方法處理該函數即可。

這里舉個例子，先隨便設一個非線性微分方程組，並給定初值：

 

然后定義最小化函數：

最后用之前介紹的非線性最優化方法解決。

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');

test1.m：

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)];

結果如下：

上面這個結果還算可以。

不過由於是非線性微分方程組，參數差一點就可能導致系統后續差別越來越大，所謂混沌與蝴蝶效應。

所以大多數擬合的結果類似下面或更糟：

注：

　　后來我又看了一下，其實這里還是屬於線性系數的，因為a,b,c,d四個系數並沒在非線性函數中。

　　非線性系數我做了幾次試驗，效果不甚理想，后面有時間再寫一篇。