ode45的案例
simulink仿真設置
一、算法設置
1.變步長(Variable—Step)求解器
可以選擇的變步長求解器有：ode45，ode23，ode113，odel5s，ode23s和discret．缺省情況下，具有狀態的系統用的是ode45；沒有狀態的系統用的是discrete。
1)ode45基于顯式Runge—Kutta(4，5)公式，Dormand—Prince對．它是—個單步求解器(solver)。也就是說它在計算y(tn)時，僅僅利用前一步的計算結果y(tn-1)．對于大多數問題．在第一次仿真時、可用ode45試一下。
2)ode23是基于顯式Runge—Kutta(2，3)．Bogackt和Shampine對．對于寬誤差容限和存在輕微剛性的系統、它比ode45更有效一些．ode23也是單步求解器。
3)odell3是變階Adams-Bashforth—Moulton PECE求解器．在誤差容限比較嚴時，它比ode45更有效．odell3是一個多步求解器，即為了計算當前的結果y(tn），不僅要知道前一步結果y(tn-1)，還要知道前幾步的結果y(tn-2)，y(tn-3)，…；
4)odel5s是基于數值微分公式(NDFs)的變階求解器．它與后向微分公式BDFs(也叫Gear方法)有聯系．但比它更有效．ode15s是一個多步求解器，如果認為一個問題是剛性的，或者在用ode45s時仿真失敗或不夠有效時，可以試試odel5s。 odel5s是基于一到五階的NDF公式的求解器．盡管公式的階數越高結果越精確，但穩定性會差一些．如果模型是剛性的，并且要求有比較好的穩定性，應將最大的階數減小到2．選擇odel5s求解器時，對話框中會顯示這一參數． 可以用ode23求解器代替。del5s，ode23是定步長、低階求解器。
展開 二分法+打靶法解微分方程
'r','s','flag','m','n','h');
s_end=5;
sp=1;
rsa=0.1;
rsb=1;
d=1;
[r1,s1]=ode45(fun,[eps,R],[rsa;0],[],m,n,h);
[r2,s2]=ode45(fun,[eps,R],[rsb;0],[],m,n,h);
figure;
subplot(121);
plot(r1,s1(:,1));hold on;
plot(r2,s2(:,1),'r');
legend('for rsa','for rsb',0);
xlim([0,R]);
tg=title(['rsa=',num2str(rsa),', rsb=',num2str(rsb)]);
subplot(122);
pa=plot(r1,s1(:,1));hold on;
pb=plot(r2,s2(:,1),'r');
pc=plot(r2,s2(:,1)*0.6,'k');
legend('for rsa','for rsb',0);
xlim([0,R]);
tg=title(['rsa=',num2str(rsa),', rsb=',num2str(rsb)]);
while abs(s1(end,1)-s2(end,1))>1e-4; % 控制精度
[r1,s1]=ode45(fun,[eps,R],[rsa;0],[],m,n,h);
[r2,s2]=ode45(fun,[eps,R],[rsb;0],[],m,n,h);
rs=[rsa+rsb]/2;
[r,s]=ode45(fun,[eps,R],[rs;0],[],m,n,h);
if s(end,1)>5; % 二分法部分
rsb=rs;
else
rsa=rs;
展開 電動汽車仿真系列-基于Simulink的防抱死制動系統
4、ABS的數學模型
防抱死制動系統的車輪轉速Simulink框圖：
整體系統框圖：
5、仿真實例
車輛初始速度：??0= 88????/??；
車輪半徑：???? = 1.25????；
車輛質量：?? = 50；
最大制動扭矩：????=1500??????*????；
液壓滯后：????=0.01 ??；
轉動慣量：????= 5????4
其中，
求解器選擇：
ode45
Solve nonstiff differential equations —medium order method
Syntax
[t,y] = ode45(odefun,tspan,y0)
[t,y] = ode45(odefun,tspan,y0,options)
[t,y,te,ye,ie] =ode45(odefun,tspan,y0,options)
sol = ode45(___)
數據檢查器選擇：
仿真結果：
模型自?。?鏈接：https://pan.baidu.com/s/17DdK3_UO1HefVRVakyh1bg
提取碼：j78t
展開 MATLAB繪制分岔圖，相平面圖以及龐加萊截面圖
,'markersize',1);end
二、MATLAB函數ode45對duffing方程進行數值求解，生成相平面圖和poincaré 截面圖
deq=@(t,x)[x(2);x(1)-0.3*x(2)-(x(1))^3+0.37*cos(1.25*t)]; options=odeset('RelTol',1e-4,'AbsTol',1e-4);
[t,xx]=ode45(deq,[0 200],[1,0],options);
plot(xx(:,1),xx(:,2),'k-','linewidth',1.5);
fsize=15;
axis([-2 2 -2 2])
xlabel('x','Fontsize',fsize);
ylabel('y','Fontsize',fsize);
title('diffing,Gamma=0.37')
clear Gamma=0.37;
deq=@(t,x)[x(2);x(1)-0.3*x(2)-(x(1))^3+Gamma*cos(1.25*t)]; options=odeset('RelTol',1e-4,'AbsTol',1e-4); [t,xx]=ode45(deq,0:(2/1.25)*pi:(4000/1.25)*pi,[1,0]); plot(xx(:,1),'.'
展開 
小球在彈簧頂端的木塊上的彈性跳動問題之代碼分析
3.代碼分析
3.1 解微分方程基本知識
微分方程分為常微分方程（ODE）與偏微分方程（PDE），簡單來說系數為常數的為ODE，而系數為變量的為PDE。目前Matlab提供了多種解微分方程的結算器，如ode113、ode15i、ode23s、ode3t等，以后有機會將系統的闡述各種算法的差異以及其各種擅長的領域，現在只需要記住，一般默認ode45算法即可，Matlab公司的建議是，如果不清楚系統的具體情況，建議先選擇ode45。常微分解法器的輸入輸出格式有：
[t,Y] = solver(odefun,tspan,y0)；
[t,Y] = solver(odefun,tspan,y0,options)；
[t,Y,TE,YE,IE] = solver(odefun,tspan,y0,options)；
上述指令的“solver”所指的就是ode45、ode45, ode23, ode113, ode15s, ode23s, ode23t, or ode23tb等，其他參數的含義介紹如下：
a) odefun 用以計算微分程右邊的函數，所有解法器處理的系統方程形式均是y’=dydt=f(t,y)的形式，通常可由odefile模板改寫而成；
b) tspan 為積分區間的向量[t0,tf],解法器設定初值在tspan(1)，然后由tspan(1)積分至tspan(end)。如果需要得到特定時間所對應的解，可以用tspan[t0,t1,…,tf]的形式。
展開 如何將vc和matlab(simulink)接口的例子
舉個例子，若為一剛性系統，即便我的程序選用了ode15s而默認為ode45,本次仿真確實用了ode15s解，仍然會報警說應該用剛性解法。不過對于剛性系統，ode45可不好使，因此從仿真效果上可以認定我們設定的剛性解法奏效了。大家不要被表面現象迷惑。
注意，要將仿真模塊放到matlab訪問的目錄下
希望能給大家帶來幫助。
246012-SimulinkVC.rar
[轉帖]繪制龐加萊截面圖
't','x','flag','betaa','F','v');
% Poincare_section[繪制龐加萊截面圖]
[t,x]=ode45(Poin,[0,2800],[0,1.5,0],[],betaa,F,v);
x(:,2)=mod(x(:,2),2*pi)-pi;
phi0=pi*2/3; % 選擇phi=2*pi/3這個截面
for k=1:round(max(x(:,3))/2/pi);
d=x(:,3)-(k-1)*2*pi-phi0;
[P,K]=sort(abs(d));
x1l=x(K(1),1);
x1r=x(K(2),1);
x2l=x(K(1),2);
x2r=x(K(2),2);
x3l=x(K(1),3);
x3r=x(K(2),3);
if abs(P(1))+abs(P(2))<3e-16;
X1(k)=x1l;
X2(k)=x2l;
else
Q=polyfit([x3l,x3r],[x1l,x1r],1);
X1(k)=polyval(Q,(k-1)*2*pi-phi0);
Q=polyfit([x3l,x3r],[x2l,x2r],1);
X2(k)=polyval(Q,(k-1)*2*pi-phi0);
end
end
plot(X1,X2,'.');
xlabel('\theta','fontsize',14);
ylabel('d\theta/dt','fontsize',14);
展開 模態疊加法和Runge-Kutta方法解動力學方程的區別
*sin(3*t)
%定義仿真時間和采樣點數
t=0:0.01:100;
%對結果進行fft變換
u1=eval(x(1));
u2=eval(x(2));
u3=eval(x(3));（省去后面的畫圖和FFt變換部分），
下面是調用ode45函數的代碼
%test4.m
function f=test4(t,y);
m1=1;m2=1;m3=2;k=1;
U=[0 1 0 0 0 0;
-3 0 1 0 0 0;
0 0 0 1 0 0;
1 0 -2 0 1 0;
0 0 0 0 0 1;
0 0 1 0 -3 0];
f=U*y+[0 sin(3*t) 0 0 0 0]';
%test4Result.m
[t,y]=ode45('test4',[0:0.01:100],[0 0 0 0 0 0]');
u1=y(:,1);
u2=y(:,3);
u3=y(:,5);
(后面省去畫圖和fft變換部分)
展開 用matlab解含分段函數的一階微分方程
'1.18*sin(10*t/pi)-u/6.7)/0.047'],'t','u');
[t,u]=ode45(fun,[0,10],[0]);
plot(t,u)
說明g這樣表示的:
gt=(sin(10*t/pi)>0)*1.18*sin(10*t/pi);
感謝蘿卜網友
解矩陣微分方程組一例
't','x','flag','n','M','K');
[t,x]=ode45(Df,[0,10],rand(n,1),[],n,M,K);
plot(t,x(:,1:n));
for k=1:n;
eval(['Le',num2str(k),'=[''X',num2str(k),'''];']);
end
ss='Le1';
for k=2:n;
ss=[ss,',Le',num2str(k)];
end
eval(['legend(',ss,',0);']);
感謝蘿卜網友
展開 [轉貼] 非線性微分方程的求解
非線性微分方程表達式的函數：
function vprime= aadlx(t,x)
% x(1)為位移 x(2)為速度 分段函數
if x(1)>1
vprime=[x(2);1.1+0.1*cos(0.5*t)+0.075*sin(0.5*t)-0.02*x(2)-(1+0.1*cos(0.5*t))*x(1)];
elseif -1<=x(1)<=1
vprime=[x(2);0.1+0.075*sin(0.5*t)-0.02*x(2)];
else
vprime=[x(2);-0.9-0.1*cos(0.5*t)+0.075*sin(0.5*t)-0.02*x(2)-(1+0.1*cos(0.5*t))*x(1)];
end
運行以下內容即可得到結果：
options=odeset('RelTol',1e-4,'AbsTol',[1e-6]);
tspan=[0,900];
[t,x]=ode45('aadlx',tspan,[0;0],options);
u1=x(:,1);
u2=x(:,2);
figure('unit','normalized','color',[1,1,1]);
h=get(gcf);
set(gcf,'Name','1','numbertitle','off');
plot(t,u1)
title('圖1')
xlabel('時間t');ylabel('位移x');
grid on
figure('unit','normalized','color',[1,1,1]);
h=get(gcf);
set(gcf,'Name','2','numbertitle','off');
plot(u1,u2)
title('圖2')
xlabel('位移x');ylabel('速度dx');
grid on
摘自 振動論壇
展開 
Python方程組獲取雅克比矩陣和海塞矩陣
前面講過當我們處理一個常微分方程組（一般對應于一個物理系統的求解）的時候，可以直接采用matlab中的ode函數比如ode15s和ode45，python的scipy中也有對應的函數像odeint和odebvp。但是，如果你想硬剛，自己去實現這些東西的話你就需要造很多的輪子，這里給出一個可以根據你輸入的代數方程組獲取雅克比矩陣和海塞矩陣的函數，并且我還特地將他做成了一個類哦，你可以直接進行實例化然后調用，非常方便，而且我還給出了測試函數。特別適合想硬剛底層算法的鐵子。
Wolfram 光學解決方案
創建光學系統的設計、曲線擬合或數據分析的互動工具，提供視覺反饋使得創新儀器的調試檢測變得容易
Code V 和 Zemax 不提供個性化的交互工具
利用完全自動的精度控制以及任意精度算法，在光學模型的計算中得出準確的結果
Matlab 和其他依賴于機械算術的系統由于缺乏數值準確度可能會出現重大錯誤
用戶可選擇所需的過程式、函數式和規則式的編程范例，使得新算法模式的建立快于其他軟件
Code V 和 Zemax 使用過程式語言
導入或獲取數據、分析數據以及遞交結果都在一個文檔中進行，無需使用多個應用程序
Wolfram 特有技術
高度優化了的超級函數分析方程，自動選擇正確的算法，以便快速得出準確結果 —— 有時為了進一步優化的需要，中途改變算法
其他計算系統要求用戶手動分析自己的方程，來確定要應用哪一個函數——例如，在 Mathematica 中您只需要使用 NDSolve 的地方，在 Matlab 中您必須要從 ode45、ode23、ode113、ode15s、bvp4c、pdepe 等中做出正確選擇，否則就會有得到錯誤結果的可能
主要功能
Wolfram技術包括用于計算、建模、可視化、開發和部署的數千種內置函數?
光學領域的專業功能
展開 matlab安裝、運行與其他問題集錦
(Type "warning off MATLAB:ode45:IntegrationTolNotMet" to suppress this warning.)
在matlab中是什么錯誤？
A： 迭代達不到所給定解的誤差。
Q21：matlab中排列組合的函數是什么？
A： nchoosek or something like that
Q22：matlab 7.0 不能在64位的cpu下運行？
A： matlab 應該是依賴于自己的虛擬機的。但是好像這個虛擬機是在IA32 里面作出來的，
所以，應該找個帶64 位的java 虛擬機替換原來的，不過不一定能行。
Q23：matlab中常用對數函數的怎么調用？
A help log
log( x ) / log( 10 ) surely works
log10( x ) ...
展開 小球在彈簧頂端的木塊上的彈性跳動問題之問題描述
PS.今天給出這個案例的全部代碼，最近將對這些代碼進行剖析，敬請期待：
%%以下是主函數，當主函數1與子函數1,2,3在同一目錄下，運行主
%%函數即可
h0=50;
k=60;
m1=20;m2=50;
tstart=0;
tfinal=1000;
y0=[h0;0;0;0];
tout=tstart;
yout=y0';
options=odeset('Events','on');
for i=1:25
[t,y,event]=ode45('xqythkfun',[tstart:0.03:tfinal],y0,options);
tout=[tout;t(2:end)];
yout=[yout;y(2:end,:)];
y0(1)=y(end,1);y0(2)=y(end,2);
v10=y(end,3);v20=y(end,4);
y0(3)=(-m2*v10+2*m2*v20+m1*v10)/(m2+m1);
y0(4)=(2*m1*v10+m2*v20-v20*m1)/(m2+m1);
tstart=t(end);
end
figure
ylabel('高度');
xlabel('時間');
hold on
plot(tout,yout(:,1),tout,yout(:,2));
legend('小球','木板');
figure
axis([-1 1 -50 h0+10]);
axis off
hold on
yt1=-45:0.3:0;
xt1=0.06*sin(yt1);
tanhuang=line(xt1,yt1,'color','k','erasemode','xor','linewidth',2);
qiu=
展開 