1樓:渠振鍾離傲冬
3.6.2
龍格-庫塔方法
改進的尤拉法比尤拉法精度高的原因在於,它在確定平均斜率時,多取了一個點的斜
率值。這樣,如果我們在[xi,x(i+1)]上多取幾個點的斜率值,然後對它們作線性組合得到平均
斜率,則有可能構造出精度更高的計算方法。這就是龍格-庫塔法的基本思想。龍格-庫塔
法可看作是尤拉法思想的提高,屬於精度較高的單步法。
龍格-庫塔法是求解常微分方程初值問題的最重要的方法之一。matlab中提供了幾
個採用龍格-庫塔法來求解常微分方程的函式,即ode23,ode45,ode113
,ode23s
,ode15s
等,其中最常用的函式是
ode23(
二三階龍格-庫塔函式)和ode45(
四五階龍格-庫塔函式),
下面分別對它們進行介紹。
1.二三階龍格-
庫塔函式(ode23)
函式ode23
的呼叫格式如下:
(1)[t,y]=ode23('f',tspan,y0)
輸入引數中的'f'
是一個字串,表示微分方程的形
式,也可以是f(x
,y)的m檔案。tspan=[t0
tfinal]表示積分割槽間,y0表示初始條件。
函式ode23
表示在初始條件
y0下從
t0到tfinal
對微分方程
'(,)
yfty
=進行積分。函式
f(t,
y)必須返回一列向量,兩個輸出引數是列向量
t與矩陣
y,其中向量
t包含估計響應
的積分點,而矩陣
y的行數與向量
t的長度相等。向量
t中的積分點不是等間距的,這是
為了保持所需的相對精度,而改變了積分演算法的步長。為了獲得在確定點t0,t1,
"的解,
tspan=[t0
t1tfinal]
。需要注意的是:tspan中的點必須是單調遞增或單調遞減的。
(2)[t,y]=ode23('f',tspan,y0,options)
其中,引數
options
為積分引數,它可由函
數odeset
來設定。options引數最常用的是相對誤差『reltol』(
預設值是
1e-3)和絕對誤差
『abstol』(預設值是
1e-6),其他引數同上。
(3)[t,y]=ode23('f',tspan,y0,options,p1,p2,…)
引數p1,p2,
…可直接輸入到函式
f中去.如
f(t,y,flag,p1,p2,…)。如果引數
options為空,則輸入
options=[
]。也可
以在ode檔案中(可參閱
odefile函式)指明引數
tspan、y0和options的值。如果參
數tspan
或y0是空,則ode23函式通過呼叫ode檔案[tspan,
y0,options]
=f([
],[],
'init
')來獲得
ode23函式沒有被提供的自變數值。如果獲得的自變數表示空,則函
數ode23會忽略,此時為
ode23('f')。
(4)[t,y,te,ye,ie]=ode23('f',tspan,y0,options)
此時要求在引數
options
中的事件屬性設為'on'
,ode檔案必須被標記,以便
p(t,y,'events')
能返回合適的資訊,詳細可參
閱函式odefile。輸出引數中的
te是一個列向量,矩陣
ye的行與列向量
te中元素相
對應,向量
ie表示解的索引。
2.四五階龍格-
庫塔函式(ode45)
函式ode45
的呼叫格式同
ode23
相同,其差別在於內部演算法不同。如果'f'
為向量函式,
則ode23
和ode45
也可用來解微分方程組。
【例3.47
】分別用二三階龍格-庫塔法和四五階龍格-庫塔法解常微分方程的初值問題:
解:先將微分方程寫成自定義函式
exam2fun.m
function
f=exam2fun
(x,y)
f=-y-x*y.^2;
f=f(:);
然後在命令視窗輸入以下語句:
>>[x1,y1]=ode23('exam2fun',[0:0.1:1],1)x1=
00.1000
0.2000
0.3000
0.4000
0.5000
0.6000
0.7000
0.8000
0.9000
1.0000y1=
1.0000
0.9006
0.8046
0.7144
0.6314
0.5563
0.4892
0.4296
0.3772
0.3312
0.2910
>>[x2,y2]=ode45('exam2fun',[0:0.1:1],1)x2=
00.1000
0.2000
0.3000
0.4000
0.5000
0.6000
0.7000
0.8000
0.9000
1.0000y2=
1.0000
0.9006
0.8046
0.7144
0.6315
0.5563
0.4892
0.4296
0.3772
0.3312
0.2910
2樓:du知道君
function dx=myfun(t,x)dx(1)=20/7*x(1)-x(2)*x(3);dx(2)=-10*x(2)+x(1)*x(3);dx(3)=-4*x(3)+x(1)*x(2);dx=dx(:); 這部分儲存為m函式檔案 命令列執行 >> x0=[3,-4,2]; t0=0.001:
0.001:20; [t,x]=ode45('myfun',[0.
001,20],x0); %ode45會自動調整步長 plot(t,x) legend('x','y','z') >> 結果
3樓:匿名使用者
matlab並沒有限制其初始值必須為0,起點可以是你定的任何數,matlab中一般以t_0表示
matlab ode45用法
4樓:大野瘦子
用法:[t,y] = ode45(odefun,tspan,y0)1、odefun 是函式
控制代碼,可以是函式檔名,匿名函式控制代碼或行內函數名。
2、tspan是區間 [t0 tf] 或者一系列散點[t0,t1,...,tf]。
3、y0是初始值向量。
4、t返回列向量的時間點。
5、y返回對應t的求解列向量。
算例程式:
function testode45
tspan=[3.9 4.0]; %求解區間y0=[8 2]; %初值
[t,x]=ode45(@odefun,tspan,y0);
plot(t,x(:,1),'-o',t,x(:,2),'-*')legend('y1','y2')
title('y'' ''=-t*y + e^t*y'' +3sin2t')
xlabel('t')
ylabel('y')
function y=odefun(t,x)y=zeros(2,1); % 列向量
y(1)=x(2);
y(2)=-t*x(1)+exp(t)*x(2)+3*sin(2*t); %常微分方程公式
endend
5樓:機智的煎餅
ode45,常微分方程的數值求解。matlab提供了求常微分方程數值解的函式。當難以求得微分方程的解析解時,可以求其數值解。matlab ode45用法如下:
[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(odefun,[t0tf],y0...)
[t,y] = ode45(odefun,tspan,y0)
odefun 是函式控制代碼,可以是函式檔名,匿名函式控制代碼或行內函數名
tspan 是區間 [t0 tf] 或者一系列散點[t0,t1,...,tf]
y0 是初始值向量
t 返回列向量的時間點
y 返回對應t的求解列向量
[t,y] = ode45(odefun,tspan,y0,options)
options 是求解引數設定,可以用odeset在計算前設定誤差,輸出引數,事件等
[t,y,te,ye,ie] =ode45(odefun,tspan,y0,options)
在設定了事件引數後的對應輸出
te 事件發生時間
ye 事件發生時之答案
ie 事件函式消失時之指標i
sol =ode45(odefun,[t0 tf],y0...)
sol 結構體輸出結果
ode的作用
ode是matlab專門用於解微分方程的功能函式。該求解器有變步長(variable-step)和定步長(fixed-step)兩種型別。
不同型別有著不同的求解器,其中ode45求解器屬於變步長的一種,採用runge-kutta演算法;其他採用相同演算法的變步長求解器還有ode23。
ode45表示採用四階-五階runge-kutta演算法,它用4階方法提供候選解,5階方法控制誤差,是一種自適應步長(變步長)的常微分方程數值解法,其整體截斷誤差為(δx)^5。
解決的是nonstiff(非剛性)常微分方程。ode45是解決數值解問題的首選方法,若長時間沒結果,應該就是剛性的,可換用ode15s試試。
6樓:ieio啊
ode45表示採用四階-五階runge-kutta演算法,它用4階方法提供候選解,5階方法控制誤差,是一種自適應步長(變步長)的常微分方程數值解法,其整體截斷誤差為(δx)^5。解決的是nonstiff(非剛性)常微分方程。
ode45語法:
[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(odefun,[t0tf],y0...)
[t,y] = ode45(odefun,tspan,y0)
odefun 是函式控制代碼,可以是函式檔名,匿名函式控制代碼或行內函數名
tspan 是區間 [t0 tf] 或者一系列散點[t0,t1,...,tf]
y0 是初始值向量
t 返回列向量的時間點
y 返回對應t的求解列向量
[t,y] = ode45(odefun,tspan,y0,options)
options 是求解引數設定,可以用odeset在計算前設定誤差,輸出引數,事件等
[t,y,te,ye,ie] =ode45(odefun,tspan,y0,options)
在設定了事件引數後的對應輸出
te 事件發生時間
ye 事件發生時之答案
ie 事件函式消失時之指標i
sol =ode45(odefun,[t0 tf],y0...)
sol 結構體輸出結果
如何在function裡使用ode45輸出值
(1) 主程式 (test.m)
邊界值為 y(1/1.5)=alpha=0 y(1)=beta=0
用 shooting method 去解二階 ode 的邊界值問題,
解 ode 使用的指令為 ode45
(2)function (funtest1.m)
解4 條first-order initial value problems
但a 的值是要從判斷解出來的值運算後,是否有大於 1 來設定
h=0.25;
m=1.2;
si=((y/x)^2-y*y'/x+(y')^2)^0.5
if si>1
a=(si.^m-1)/(h*si)
elseif si<=1
a=0end
MATLAB關於一段ode45的程式報錯
劉復民 我給你執行了一下沒錯啊 這是一部分資料 不過你在第一個函式儲存時函式名應該是 hj把 266.7452 266.7543 266.7634 266.7725 266.7816 266.7907 266.7997 266.8088 266.8179 266.8270 266.8361 266....
關於用matlab處理影象的問題,如何將下進行邊緣提取,得到比較清晰的輪廓,最好是能尋找到偏心孔的座標
gimg rgb2gray img img是彩色影象,gimg是灰度影象 eimg edge gimg,sobel sobel邊緣檢測,eimg是二值影象 剩下的可以用數學形態 算,計算每個連通域的面積和周長,圓的周長面積滿足一定比例,滿足此比例的就是偏心孔位置 敬詩蘭 用matlab處理 除了提到...
關於matlab的quad函式,matlab中quad命令點除點乘,是什麼意思? quad x 4 1 x 2 ,
quadl q.p sin p 1,2,2 這個函式的變數是q而不是p,符號計算時,matlab自動呼叫函式確定變數,預設為最靠近x的作為變數 盧憐陽衡芬 y x 4 3 sqrt 1 2 3 2.sin x 2 true quad y,0,2 pi 按上面語句試試 x 表示函式以x為自變數,mat...