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

比如我們已經有了微分方程模型和相關數據，如何求模型的參數。

這里以SEIR模型為例子，SEIR模型可以參考之前的文章。

一般的線性方程我們可以用最小二乘來解，一般的非線性方程我們可以用LM來解。

這里是線性微分方程組，所以我們采用最小二乘來解。

關鍵是構造出最小二乘形式，微分可以通過前后數據差分的方法來求。

不過這里還有一個技巧就是如果數據前后幀間隔過大，可以先插值，再對插值后的數據差分。

如果實際測量數據抖動過大導致插值后差分明顯不能反映實際情況，可以先對數據平滑（擬合或是平均）再求差分。

先看SEIR微分方程：

寫成矩陣形式：

 

到這里就能用最小二乘來求解了。 

matlab代碼如下：

main.m：

clear all;close all;clc;

%%SEIR模型
A = [0.5 0.1 0.05 0.02];
[t,h] = ode45(@(t,x)SEIR(t,x,A),[0 300],[0.01 0.98 0.01 0]);  %[初始感染人口占比 初始健康人口占比 初始潛伏人口占比 初始治愈人口占比]
plot(t,h(:,1),'r');
hold on;
plot(t,h(:,2),'b');
plot(t,h(:,3),'m');
plot(t,h(:,4),'g');
legend('感染人口占比I','健康人口占比S','潛伏人口占比E','治愈人口占比R');
title('SEIR模型')

data=[t h];
data = data(1:3:80,:);      %間隔取一部分數據用來擬合
figure;
plot(data(:,1),data(:,2),'ro');
hold on;
plot(data(:,1),data(:,3),'bo');
plot(data(:,1),data(:,4),'mo');
plot(data(:,1),data(:,5),'go');

T=min(data(:,1)):0.1:max(data(:,1));        %插值處理，如果數據多，也可以不插值
I=spline(data(:,1),data(:,2),T)';
S=spline(data(:,1),data(:,3),T)';
E=spline(data(:,1),data(:,4),T)';
R=spline(data(:,1),data(:,5),T)';

plot(T,I,'r.');plot(T,S,'b.');
plot(T,E,'m.');plot(T,R,'g.');

%求微分，如果數據幀間導數變化太大，可以先平均或者擬合估計一個導數
%因為前面T是以0.1為步長，這里乘以10
dI = diff(I)*10; dI=[dI;dI(end)];       
dS = diff(S)*10; dS=[dS;dS(end)];
dE = diff(E)*10; dE=[dE;dE(end)];
dR = diff(R)*10; dR=[dR;dR(end)];

X = [zeros(length(I),1) -I.*S zeros(length(I),2);   %構造線性最小二乘方程組形式
     -E I.*S -E zeros(length(I),1);
     E zeros(length(I),2) -I;
     zeros(length(I),2) E I];
Y = [dS;dE;dI;dR];

A = inv(X'*X)*X'*Y

%用估計參數代入模型
[t,h] = ode45(@(t,x)SEIR(t,x,A),[0 300],[I(1) S(1) E(1) R(1)]);  %[初始感染人口占比 初始健康人口占比 初始潛伏人口占比 初始治愈人口占比]
plot(t,h(:,1),'r');
hold on;
plot(t,h(:,2),'b');
plot(t,h(:,3),'m');
plot(t,h(:,4),'g');

SEIR.m：

function dy=SEIR(t,x,A)
alpha = A(1);        %潛伏期轉陽率
beta = A(2);         %感染率
gamma1 = A(3);      %潛伏期治愈率
gamma2 = A(4);      %患者治愈率
dy=[alpha*x(3) - gamma2*x(1);
    -beta*x(1)*x(2);
    beta*x(1)*x(2) - (alpha+gamma1)*x(3);
    gamma1*x(3)+gamma2*x(1)];

結果：

原始參數[0.5 0.1 0.05 0.02]與模型：

擬合參數[0.499921929359668 0.100099242849855 0.0505821757746970 0.0199739921888752]與模型：