微分方程式の解を複数​の入力値があるときで​も適応させたい。

まず、作成したいと思っているプログラムの趣旨について説明いたします。

本プログラムでは、微分方程式によって得た値に応じてプロットを行います。

微分方程式で求める解は、x、y、θの3変数になります。

また、微分方程式を求めるために使用する関数はOde45です。

ode45には初期値が必要であると思いますが、ループ1回目の初期値は[0,0,0]'とし、

2回目以降の初期値は求めた3つの変数の最後の配列に格納された値を使用します。

以下に、現在作成中のプログラムを転記いたします。

clc;

clear;

close all;

%%

%パラメータ設定

P = struct

%機体中心から車輪中心までの長さ

P.l = 340;%[mm]

%タイヤ半径

P.r = 125;%[mm]

%変数

%左右車輪の回転数

r_R = [10,20];%[rpm]

r_L = [50,50];%[rpm]

%積分範囲

plotnumber = 11;

tspan = linspace(0,1,plotnumber); % 0 ~ 1までの範囲

for ii = 1:length(r_R)

%左右車輪の角速度

w_R = (r_R(ii)*2*pi)/60;%[rad/s]

w_L = (r_L(ii)*2*pi)/60;%[rad/s]

%角速度から車輪速度をv=rwで求める

v_R = (w_R * P.r);%[mm/s]

v_L = (w_L * P.r);%[mm/s]

% w,V,pの計算部分

% wは旋回中心から見たときの角速度

w = (v_R - v_L)/(2*P.l) %[rad/s]

% Vは本体の速度

V = (v_R + v_L)/2 %[mm/s]

% pは機体中心から見た旋回中心までの距離

if w_R > w_L

p = (P.l*(v_R + v_L))/(v_R - v_L) %[mm]

%左車輪中心

p_lc = p - P.l

%右車輪中心

p_rc = p + P.l

elseif w_R < w_L

p = (P.l*(v_R + v_L))/(v_R - v_L) %[mm]

%左車輪中心

p_lc = p + P.l

%右車輪中心

p_rc = p - P.l

end

%%

%初期値 (x,y) = (0,0),進行方向θ=0

if ii == 1

q0 = [0,0,0]';

else

q0 = [q(plotnumber,1),q(plotnumber,2),q(plotnumber,3)]';

end

%%

%微分方程式を解く

%ode45は(関数がtで変化することを記載 関数,tの範囲,初期値)で成り立つ

[t,q] = ode45(@(t,q) f(t,q,P,w_R(ii),w_L(ii)),tspan,q0);

%%

%描画

figure

title('two wheel robot center and orbit by ODE45');

hold on

grid on

for j = 1:length(t)

%機体中心を描画

plot(q(j,1) , q(j,2),'ro')

end

for j = 2:length(t)

%機体中心の差分を描画

plot([q(j,1) , q(j-1,1)] , [q(j,2) , q(j-1,2)],'Color','red','LineStyle','-')

end

axis equal

end

hold off

%%

function dydt = f(t,q,P,w_R,w_L)

%パラメータ設定

r = P.r;

l = P.l;

%状態方程式

dydt = [r/2*(w_R + w_L)*cos(q(3))

r/2*(w_R + w_L)*sin(q(3))

r/(2*l) * (w_R - w_L) ];

end

また、このプログラムを実行した際のエラーについても同様に転記いたします。

インデックスが配列要素数 (1) を超えています。

エラー: two_wheel_robot_orbit>@(t,q)f(t,q,P,w_R(ii),w_L(ii)) (行 57)

[t,q] = ode45(@(t,q) f(t,q,P,w_R(ii),w_L(ii)),tspan,q0);

エラー: odearguments (行 90)

f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.

エラー: ode45 (行 115)

odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

エラー: two_wheel_robot_orbit (行 57)

[t,q] = ode45(@(t,q) f(t,q,P,w_R(ii),w_L(ii)),tspan,q0);

よろしくお願いいたします。

Answer 1

[Atsushi Ueno]

発生しているエラーは下記変更で解消します。

%  [t,q] = ode45(@(t,q) f(t,q,P,w_R(ii),w_L(ii)),tspan,q0); % 変更前
%  [t,q] = ode45(@(t,q) f(t,q,P,w_R,w_L),tspan,q0); % 変更後

r_Rは2要素のベクトルであるのに対し、r_Rの1要素を選択して下記の通り演算する

%  w_R = (r_R(ii)*2*pi)/60;%[rad/s]

w_Rは1要素のスカラ値です。その上で下記命令を実行すると、

%  [t,q] = ode45(@(t,q) f(t,q,P,w_R(ii),w_L(ii)),tspan,q0);

ii=1の時はw_R(1)にアクセス出来ますが、ii=2の時はw_R(2)が存在しないのでエラーになります。

(追記)＞1つの図に計算結果を繰り返し出力するためにはどのように行えばよろしいでしょうか。

figureとhold onをforループの外に出せば同じプロット画面に繰り返し出力するようになります。

clc; clear; close all;
P = struct;
P.l = 340;%[mm]%機体中心から車輪中心までの長さ
P.r = 125;%[mm]%タイヤ半径
r_R = [10,20];%[rpm]%左右車輪の回転数
r_L = [50,50];%[rpm]%左右車輪の回転数
plotnumber = 11;%積分範囲
tspan = linspace(0,1,plotnumber); % 0 ~ 1までの範囲
figure
title('two wheel robot center and orbit by ODE45');
hold on
grid on
for ii = 1:length(r_R)
    w_R = (r_R(ii)*2*pi)/60;%[rad/s] %左右車輪の角速度
    w_L = (r_L(ii)*2*pi)/60;%[rad/s]
    v_R = (w_R * P.r);%[mm/s] %角速度から車輪速度をv=rwで求める
    v_L = (w_L * P.r);%[mm/s]
    w = (v_R - v_L)/(2*P.l); %[rad/s]% wは旋回中心から見たときの角速度
    V = (v_R + v_L)/2; %[mm/s]% Vは本体の速度
    if w_R > w_L % pは機体中心から見た旋回中心までの距離
        p = (P.l*(v_R + v_L))/(v_R - v_L); %[mm]
        p_lc = p - P.l; %左車輪中心
        p_rc = p + P.l; %右車輪中心
    elseif w_R < w_L
        p = (P.l*(v_R + v_L))/(v_R - v_L); %[mm]
        p_lc = p + P.l; %左車輪中心
        p_rc = p - P.l; %右車輪中心
    end
    if ii == 1
        q0 = [0,0,0]';%初期値 (x,y) = (0,0),進行方向θ=0
    else
        q0 = [q(plotnumber,1),q(plotnumber,2),q(plotnumber,3)]';
    end
    %微分方程式を解く
    %ode45は(関数がtで変化することを記載 関数,tの範囲,初期値)で成り立つ
    [t,q] = ode45(@(t,q) f(t,q,P,w_R,w_L),tspan,q0);
    %描画
    for j = 1:length(t)
        plot(q(j,1) , q(j,2),'ro')%機体中心を描画
    end
    for j = 2:length(t)
        plot([q(j,1) , q(j-1,1)] , [q(j,2) , q(j-1,2)],'Color','red','LineStyle','-')%機体中心の差分を描画
    end
    axis equal
end
hold off

function dydt = f(t,q,P,w_R,w_L)
%パラメータ設定
r = P.r;
l = P.l;
%状態方程式
dydt = [r/2*(w_R + w_L)*cos(q(3))
        r/2*(w_R + w_L)*sin(q(3))
        r/(2*l) * (w_R - w_L)       ];
end

[Naoya Tanabe]

ご回答ありがとうございます。

ode45内のw_Rとw_Lをご指摘の通りに変更したところ、エラーが出ずに最後まで出力することができました。

しかし、このプログラムを実行すると、プロット図が複数出てきてしまいます。

(ループ1回目のfigureと2回目のfigureが出力される)

1つの図に計算結果を繰り返し出力するためにはどのように行えばよろしいでしょうか。

[Atsushi Ueno]

hold onです

新しいプロットを追加するときに現在のプロットを保持 - MATLAB hold - MathWorks 日本

[Atsushi Ueno]

質問文のプログラムを変更して、1つの図に計算結果を繰り返し出力するための方法を回答に追記しました。

[Naoya Tanabe]

追記していただいたプログラムを実行することで、希望していたグラフをプロットすることができました！

大変ありがとうございました！