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

Question

まず、作成したいと思っているプログラムの趣旨について説明いたします。
本プログラムでは、微分方程式によって得た値に応じてプロットを行います。
微分方程式で求める解は、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);

よろしくお願いいたします。

Atsushi Ueno · Accepted Answer

発生しているエラーは下記変更で解消します。
%  [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

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

4 个评论
显示 2更早的评论隐藏 2更早的评论

更多回答（0 个）

类别

产品

版本

标签

Community Treasure Hunt

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

0 个评论 显示 -2更早的评论 隐藏 -2更早的评论

采纳的回答

4 个评论 显示 2更早的评论 隐藏 2更早的评论

更多回答（0 个）

类别

产品

版本

标签

另请参阅

Community Treasure Hunt

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

4 个评论
显示 2更早的评论隐藏 2更早的评论