ode求解器的事件 Event 屬性

2021-07-03 15:56:52 字數 3105 閱讀 4556

檢測事件

matlab微分方程如何設定變數的範圍 

如dy1=y2

dy2=y1+1

其中y1的範圍為0

【解】

m檔案:

function [value,isterminal,direction] = events1(t,y)

value = y(1)-4;

isterminal= 1;

direction = 0;

命令視窗:

dy = @(t,y) [y(2);y(1) + 1];

options=odeset('events',@events1);

[t,y] = ode45(dy,[0 12],[0 1],options);

plot(t,y)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%function [value,isterminal,direction]=events(t,x)

% 事件檢查函式,此時需要做的是過零點檢測

% ode45函式自動檢查當value=0是否成立

% 如果我們要求檢測y=0的點,設定value=y

% 這裡我們要檢測y=4,那麼就設定value=y-4

% isterminal

檢測到指定條件時,是否終止ode45函式的執行

% 1表示終止,0表示繼續

% 在我們這個問題上,我們只要檢測到零點時就停止程式

% direction:

value過零點檢測的方向

% -1表示由正到負,+1表示由負到正

% 對於我們這個問題,當然是由正到負

%.........................

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

乙個例子和說明:

這是乙個用到簡單微分方程的物理情景

乙個質量m=100kg的物體從高處豎直落下,加速度會受到空氣阻力的影響,這裡簡單的認為重力加速度g=9.8不變,空氣阻力f=k*v^2 ,簡單起見,k=1。只考慮豎直方向速度v,且速度,加速度,豎直位移都以向下為正。初速度,初位移都為0;那麼有以下微分方程:

dy/dt=v

dv/dt=9.8-1*v^2/m

m=100,v0=y0=0

然後我用matlab的ode45函式求這個微分方程的數值解

先編寫函式

function dx=fun(t,x)

% x(1)表示下落的距離y(向下為正),x(2)表示下落速度v(向下為正)

k=1; % k=1為表示空氣阻力的乙個常量,這裡簡化空氣阻力f=k*v^2

m=100;% m為質量=100kg

dx=zeros(2,1);

dx(1)=x(2); % 下落距離對時間的導數=速度

a=9.8-(k*x(2)^2)/m;% a加速度(向下為正)=重力加速度 - 空氣阻力產生的加速度

dx(2)=a; % 速度對時間的導數=加速度

end現在我想要得到t=15s時的位移和速度

那麼輸入

[t,x]=ode45('fun',[0,15],[0 0]);

返回的x中的最後一列就是我想要的值;

但假如我想知道當豎直向下的位移剛好=100公尺時的時間和速度,那該怎麼辦?現在我的做法是先將解乙個充分大的時間,然後在裡面找位移在100兩側的時間和速度,再通過插值得到位移剛好=100時的時間和速度。但這樣很麻煩,也不見得準確,matlab有什麼自帶的語句能實現這個功能嗎?或是有什麼更好的方法?

在不知道結果時間的時候是需要先設定乙個比較大的時間範圍計算的

但是並不需要將整個範圍的結果都算出來再插值

這個時候可以設定觸發事件函式在一定條件下停止計算

用odeset可以為ode45求解器設定觸發事件的函式

詳細的用法要仔細檢視matlab的幫助檔案,這裡我以你的題目,舉個例子

微分方程還是用你的函式fun

在用ode45解方程之前,再寫乙個函式:事件觸發函式eventfun,

它的格式固定要返回三個值,這三個值的意思是

當第乙個值vaule的值到達0時,時間會觸發

而根據第二個值isterminal設定,觸發時間會否終止求解

第三個值是設定觸發過0的方向

function [value,isterminal,direction] = eventfun(t,x)

value=x(1)-100; %觸發時間,當其值為0的時候,時間會觸發

isterminal=1; %設為1時會,觸發時間會停止求解器,設0時觸發不影響工作

direction=1; %觸發方向設1時是上公升觸發,設-1是下降觸發,設0是雙向觸發

end寫好fun和eventfun之後,你就可以呼叫ode45解方程了

但用ode45之前記得先用odeset,將觸發函式加入哦

op=odeset('events',@eventfun);

[t,x,tend,xend,evennum]=ode45(@fun,[0,15],[0 0],op);

這樣你看看,到達100公尺時,求解器就停住了

細心的你注意,ode45多返回了tend,xend,evennum三個引數

第乙個tend是觸發事件發生的時間

第二個xend是觸發時間發生時刻的x

第三個evenum是標識觸發事件的編號

由於這裡只設定了乙個觸發事件,所以編號肯定是1

其實你也可以將該時間isterminal設為0

那麼求解器不會因為觸發事件而停止15秒內的解都會算出

你仍舊可以根據tend 和 xend得到到達100公尺時的時間和狀態

額外提示,有時候你不知道到底取多大的時間範圍才能夠等到你要的觸發時間

那麼你可以用迴圈判斷的方法,先設乙個時間範圍,然後求解

到最後都沒有等到觸發事件(tend和xend都是空矩陣)

那麼適當延長求解時間區間,將上次的最後時刻和狀態作為初始狀態

再一次求解時間範圍更大的解,如此直到找到觸發事件才停止

引用:1.

2.

ode求解器的事件 Event 屬性

檢測事件 matlab微分方程如何設定變數的範圍 如dy1 y2 dy2 y1 1 其中y1的範圍為0 解 m檔案 function value,isterminal,direction events1 t,y value y 1 4 isterminal 1 direction 0 命令視窗 dy...

JavaScript事件event物件屬性

ie和ff獲取事件的不同 var e window.event e.target event 物件只在事件發生的過程中才有效。阻止事件預設行為 ie window.event.returnvalue true ff e.preventdefault 阻止事件冒泡行為 ie window.event....

event 事件的使用

event用於兩個執行緒間的協作,比如乙個執行緒得到了資料發訊號給另乙個執行緒讓它來處理 多執行緒裡的event from threading import thread,event import time,random event event def light print light is li...