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

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

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

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

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

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

不過這裡還有乙個技巧就是如果資料前後幀間隔過大，可以先插值，再對插值後的資料差分。

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

先看seir微分方程：

寫成矩陣形式：

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

matlab**如下：

main.m：

clear all;close all;clc;

%%seir模型
a = [0.5
0.10.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');

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.499921929359668 0.100099242849855 0.0505821757746970 0.0199739921888752]與模型：

matlab練習程式（常微分方程組向量場）

過去有畫過常微分方程的向量場，通過向量場能夠很形象的看出方程解的狀態。這裡用matlab也實現一下，同時對三維情況也做了乙個實現。繪製的方法就是計算方程在二維或三維某個點的方向，然後把方向歸一化，畫出歸一化的向量即可。二維微分方程組如下 三維微分方程組如下 matlab 如下 二維情況 clear ...

matlab練習程式（高階常微分方程組數值解）

這裡以三元二次常微分方程組做乙個例子，更多元更高次的都類似。比如下列方程組 x x x y z y y y x z z z x matlab 如下 main.m clear all close all clc t,x ode45 testfun 0 0.1 5 1 1 1 10.5 0 0 0.1 ...

MATLAB學習筆記 常微分方程的數值解

常微分方程數值求解的命令 求常微分方程的數值解，matlab的命令格式為 t,y solver odefun tspan,y0,options 其中solver選擇ode45等函式名，odefun為根據待解方程或方程組編寫的m檔名，tspan為自變數的區間 t0,tf 即準備在那個區間上求解，y0表...