在不知道結果時間的時候是需要先設定一個比較大的時間範圍計算的
但是並不需要將整個範圍的結果都算出來再插值
這個時候可以設定觸發事件函式在一定條件下停止計算
用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都是空矩陣)
那麼適當延長求解時間區間,將上次的最後時刻和狀態作為初始狀態
再一次求解時間範圍更大的解,如此直到找到觸發事件才停止
在不知道結果時間的時候是需要先設定一個比較大的時間範圍計算的
但是並不需要將整個範圍的結果都算出來再插值
這個時候可以設定觸發事件函式在一定條件下停止計算
用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都是空矩陣)
那麼適當延長求解時間區間,將上次的最後時刻和狀態作為初始狀態
再一次求解時間範圍更大的解,如此直到找到觸發事件才停止