functionで関数を作る方法

下記のようなやり方で

蔵本モデルの微分方程式を解いたのですが、一番最初のfunctionの所でエラーが出てきてしまいます。

function dydt=kuramotomodel(t,y)

↑

エラー: このコンテキストでは関数定義はサポートされません。コードファイルで関数を作成してください。

コードファイルで関数を作成とはどのようなやり方なのでしょうか？

やり方分かるいましたら解決策を伺いたいです。

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

function dydt=kuramotomodel(t,y) K=10; N=3; omega=0.5; domega=zeros(1,N); for k=1:length(domega) domega(k)=omega; for m=1:length(domega) domega(k)=domega(k)+K/N*sin(y(m)-y(k)); end end dydt = [domega(1); domega(2); domega(3)];

function dydt=kuramotomodel(t,y)

K=10;

N=2;

omega=0.5;

domega=zeros(1,N);

for k=1:length(domega)

domega(k)=omega;

for m=1:length(domega)

domega(k)=domega(k)+K/N*sin(y(m)-y(k));

end

end

dydt = [domega(1); domega(2)];function dydt=kuramotomodel(t,y)

Answer 1

[Kojiro Saito]

コマンドウィンドウで上記の関数を定義されようとしていませんでしょうか？

コードファイルというのはスクリプトかライブスクリプトのことで、この中で関数を書く必要があります。

MATLABのホームタブから「新規スクリプト」をクリックしてエディターを開き、

上記の関数の定義をコピペし、kuramotomodel.mの名前で保存します。

するとkuramotomodelの関数が呼び出せるようになります。

または、MATLABのホームタブから「レイアウト」をクリックしてコマンド履歴を表示させ、コマンドウィンドウでじっこうしたコマンドを選択して右クリックでスクリプトの作成をクリックしてスクリプトにコピーさせることもできます。

[Kojiro Saito]

szさんのコメントをこちらに移動させていただきます。

sz write on 13:30, 2020/07/13.

---

教えていただきありがとうございます

function でのエラーがなくなりました。

しかしdomega(k)=domega(k)+K/N*sin(y(m)-y(k));

のy が認識されないとなってしまうのですが、

蔵本モデルの微分方程式を解く解いているのですが、間違えている部分分かる方いましたら

教えていただきたいです。

何度も申し訳ありませんが、よろしくお願いいたします。

function dydt=kuramotomodel(t,y)

K=10;

N=2;

omega=0.5;

domega=zeros(1,N);

for k=1:length(domega)

domega(k)=omega;

for m=1:length(domega)

domega(k)=domega(k)+K/N*sin(y(m)-y(k));

end

end

dydt = [domega(1); domega(2)];

[t,y] = ode45(@kuramotomodel,[2/3*pi;0;4/3*pi]);

v = VideoWriter('kuramotomodel.avi');

open(v);

figure(1)

for i=1:length(y(:,1))

xlim([-1.5 1.5])

ylim([-1.5 1.5])

hold on

plot(cos(y(i,1)),sin(y(i,1)),'b.',cos(y(i,2)),sin(y(i,2)),'r.',cos(y(i,3)),sin(y(i,3)),'g.')

drawnow limitrate

frame=getframe(gcf);

writeVideo(v,frame);

clf

end

close(v);

---

[Kojiro Saito]

szさん

変数yは関数kuramotomodelの2つめの入力引数として与えるようになっています。

関数をkuramotomodelをkuramotomodel.mに書いた後、

コマンドウィンドウから呼び出すときに、変数tとyを定義してから下記のように呼び出ししていますか？

dydt=kuramotomodel(t,y)

[sz]

新規スクリプトを開いてエディタに全てコピペして、実行を押すと保存する画面が出てファイル名が最初からkuramotomodel.mとなっているのでそのまま保存したらその下部分も実行されたのですがyの部分でエラーが出てしまいました。

もし可能でしたら、コマンドウィンドウから呼び出す方法と変数t とy の定義の仕方伺いたいです。

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

[Kojiro Saito]

エディターで実行を押すと、その関数が入力引数無しで実行されますので、tとyが与えられないまま処理が進み、yの部分でエラーが出てしまいます。

こちらのドキュメントのように、実行ボタンを押し下げて「コードを入力して実行」と押して、

t=0:0.1:1; y=[2/3*pi; 0; 4/3*pi]; kuramotomodel(t, y)

のように入力して実行すると入力引数を入れて実行できます。

https://jp.mathworks.com/help/matlab/matlab_prog/run-functions-in-the-editor.html

コマンドウィンドウでも同様に

t=0:0.1:1;
y=[2/3*pi; 0; 4/3*pi];
kuramotomodel(t, y)

などで実行できます。

上記のode45の式ですが、

 [t,y] = ode45(@kuramotomodel,[2/3*pi;0;4/3*pi]);

ode45のドキュメントによると ode45(odefun, tspan, y0)のように、積分区間tspanと初期条件y0の2変数が要になりますが、[2/3*pi;0;4/3*pi]がy0に該当しますでしょうか？積分区間tspanは分かりますでしょうか？

また、蔵本モデルですが、MATLABを使ったサンプルを見付けました。

Kuramoto Model of Synchronized Oscillators (ブログ記事)

ブログ記事で紹介された蔵本モデルのソースコード

ご参考になれば幸いです。

[sz]

ありがとうございます。

下記のようにコマンドウィンドウに打ったらこのようなエラーが出てきてしまいした。

>> t=0:0.1:1;

y=[2/3*pi; 0; 4/3*pi];

kuramotomodel(t, y);

エラー: odearguments (line 21)

ode45 への 1 番目の引数が関数ハンドルの場合、tspan 引数は少なくとも 2 つの要素をもたなければなりません。

エラー: ode45 (line 115)

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

エラー: kuramotomodel (line 17)

[t,y] = ode45(@kuramotomodel,t,y);

エラー: kuramotomodel (line 19)

[t,y] = ode45(@kuramotomodel,t,y);

このようなエラーがでてから、[t y]=ode45(@kuramotomodel,t,y)とはコマンドウィンドウには書いていないのですがodeargument のような感じで書かれた新しいエディターが開きました。

自分は普段エディターにスクリプトを書いて実行を押しているのですが

そこのファイルとkuramotomodel関数は別ということなのでしょうか？

蔵本モデル関数を呼び出すにはエディター及びコマンドウィンドウにどのように書けば良いのか再度伺っても宜しいでしょうか。

理解が遅く申し訳ありません。

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

Answer 2

[sz]

教えていただきありがとうございます

function でのエラーがなくなりました。

しかしdomega(k)=domega(k)+K/N*sin(y(m)-y(k));

のy が認識されないとなってしまうのですが、

蔵本モデルの微分方程式を解く解いているのですが、間違えている部分分かる方いましたら

教えていただきたいです。

何度も申し訳ありませんが、よろしくお願いいたします。

function dydt=kuramotomodel(t,y)

K=10;

N=2;

omega=0.5;

domega=zeros(1,N);

for k=1:length(domega)

domega(k)=omega;

for m=1:length(domega)

domega(k)=domega(k)+K/N*sin(y(m)-y(k));

end

end

dydt = [domega(1); domega(2)];

[t,y] = ode45(@kuramotomodel,[2/3*pi;0;4/3*pi]);

v = VideoWriter('kuramotomodel.avi');

open(v);

figure(1)

for i=1:length(y(:,1))

xlim([-1.5 1.5])

ylim([-1.5 1.5])

hold on

plot(cos(y(i,1)),sin(y(i,1)),'b.',cos(y(i,2)),sin(y(i,2)),'r.',cos(y(i,3)),sin(y(i,3)),'g.')

drawnow limitrate

frame=getframe(gcf);

writeVideo(v,frame);

clf

end

close(v);