洛侖茲方程組：定義和數值解

洛侖茲方程組是常微分方程組(請參閱洛侖茲方程組)。對于實數常量 σ,ρ,β，方程組為

dxdt=σ(y-x)dydt=x(ρ-z)-ydzdt=xy-βz.

一個敏感系統的洛侖茲參數值是 σ=10,β=8/3,ρ=28。從 [x(0),y(0),z(0)] = [10,20,10] 開始啟動系統，查看系統從時間 0 到 100 的演變。

sigma = 10;

beta = 8/3;

rho = 28;

f = @(t,a) [-sigma*a(1) + sigma*a(2); rho*a(1) - a(2) - a(1)*a(3); -beta*a(3) + a(1)*a(2)];

xt0 = [10,20,10];

[tspan,a] = ode45(f,[0 100],xt0); % Runge-Kutta 4th/5th order ODE solver

figure

plot3(a(:,1),a(:,2),a(:,3))

view([-10.0 -2.0])

演變相當復雜。但在很短的時間區間內，它看起來有點像勻速圓周運動。繪制 [0,1/10] 時間區間內的解。

[tspan,a] = ode45(f,[0 1/10],xt0); % Runge-Kutta 4th/5th order ODE solver

figure

plot3(a(:,1),a(:,2),a(:,3))

view([-30 -70])

對 ODE 解進行圓周路徑擬合

圓周軌跡的方程組有幾個參數：

路徑與 x-y 平面間的夾角 θ(1)

平面與 x 軸向傾角間的夾角 θ(2)

半徑 R

速度 V

相對于時間 0 的時移 t0

空間 delta 中的三維位移

根據這些參數，可確定時間 xdata 的圓周軌跡的位置。

type fitlorenzfn

function f = fitlorenzfn(x,xdata)

theta = x(1:2);

R = x(3);

V = x(4);

t0 = x(5);

delta = x(6:8);

f = zeros(length(xdata),3);

f(:,3) = R*sin(theta(1))*sin(V*(xdata - t0)) + delta(3);

f(:,1) = R*cos(V*(xdata - t0))*cos(theta(2)) ...

- R*sin(V*(xdata - t0))*cos(theta(1))*sin(theta(2)) + delta(1);

f(:,2) = R*sin(V*(xdata - t0))*cos(theta(1))*cos(theta(2)) ...

- R*cos(V*(xdata - t0))*sin(theta(2)) + delta(2);

要求得洛侖茲方程組在 ODE 解的時間點處的最佳擬合圓周軌跡，請使用 lsqcurvefit。為了將參數保持在合理的范圍內，可對各參數設置邊界。

lb = [-pi/2,-pi,5,-15,-pi,-40,-40,-40];

ub = [pi/2,pi,60,15,pi,40,40,40];

theta0 = [0;0];

R0 = 20;

V0 = 1;

t0 = 0;

delta0 = zeros(3,1);

x0 = [theta0;R0;V0;t0;delta0];

[xbest,resnorm,residual] = lsqcurvefit(@fitlorenzfn,x0,tspan,a,lb,ub);

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to

its initial value is less than the value of the function tolerance.

繪制 ODE 解的時間點處的最佳擬合圓周點以及洛侖茲方程組的解。

soln = a + residual;

hold on

plot3(soln(:,1),soln(:,2),soln(:,3),'r')

legend('ODE Solution','Circular Arc')

hold off

figure

plot3(a(:,1),a(:,2),a(:,3),'b.','MarkerSize',10)

hold on

plot3(soln(:,1),soln(:,2),soln(:,3),'rx','MarkerSize',10)

legend('ODE Solution','Circular Arc')

hold off

對圓弧進行 ODE 擬合

現在修改參數 σ,β,andρ，使其與圓弧最佳擬合。為了更好地擬合，也允許初始點 [10,20,10] 發生改變。

為此，編寫函數文件 paramfun，該文件采用 ODE 擬合的參數，并計算在時間 t 上的軌跡。

type paramfun

function pos = paramfun(x,tspan)

sigma = x(1);

beta = x(2);

rho = x(3);

xt0 = x(4:6);

f = @(t,a) [-sigma*a(1) + sigma*a(2); rho*a(1) - a(2) - a(1)*a(3); -beta*a(3) + a(1)*a(2)];

[~,pos] = ode45(f,tspan,xt0);

要找到最佳參數，可使用 lsqcurvefit 來最小化新計算的 ODE 軌跡和圓弧 soln 之差。

xt0 = zeros(1,6);

xt0(1) = sigma;

xt0(2) = beta;

xt0(3) = rho;

xt0(4:6) = soln(1,:);

[pbest,presnorm,presidual,exitflag,output] = lsqcurvefit(@paramfun,xt0,tspan,soln);

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to

its initial value is less than the value of the function tolerance.

確定這種優化對參數的改變程度。

fprintf('New parameters: %f, %f, %f',pbest(1:3))

New parameters: 9.132446, 2.854998, 27.937986

fprintf('Original parameters: %f, %f, %f',[sigma,beta,rho])

Original parameters: 10.000000, 2.666667, 28.000000

參數 sigma 和 beta 大約改變 10%。

繪制修正后的解。

figure

hold on

odesl = presidual + soln;

plot3(odesl(:,1),odesl(:,2),odesl(:,3),'b')

plot3(soln(:,1),soln(:,2),soln(:,3),'r')

legend('ODE Solution','Circular Arc')

view([-30 -70])

hold off

ODE 擬合中的問題

如優化仿真或常微分方程中所述，優化器可能會因為數值 ODE 解中的固有噪聲而遇到問題。如果您懷疑解不理想(例如，退出消息或退出標志指示解可能不準確)，則可嘗試更改有限差分運算。在本示例中，我們使用了更大的有限差分步長和中心有限差分。

options = optimoptions('lsqcurvefit','FiniteDifferenceStepSize',1e-4,...

'FiniteDifferenceType','central');

[pbest2,presnorm2,presidual2,exitflag2,output2] = ...

lsqcurvefit(@paramfun,xt0,tspan,soln,[],[],options);

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to

its initial value is less than the value of the function tolerance.

在這一情況中，使用這些有限差分選項不能改進解。

disp([presnorm,presnorm2])

20.0637 20.0637

總結

以上是生活随笔為你收集整理的matlab 微分方程组参数拟合,拟合常微分方程 (ODE)的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。