combvecによる​メモリ不足および積分​計算に関しまして

Xにより7つのパラメータを遷移させた上で，解としてk3を求めるプログラムを組んだのですが，エラーが発生してしまいました．

操作としては，Xで全てのパラメータのパターンを網羅したベクトルをつくり，別のmファイル(trap2.m)から関数を呼び出し各パラメータを代入し，方程式eqnを満たすk3を求めるといったものになります．

また，方程式eqnはパラメータk1，k2，t1，t2，t3，p1，p2を与えた上で，trap2.m(台形を組み合わせた関数)の面積が０となるようなk3を求めるといったものです．

trap2.mの概形は写真にて添付させていただきました．

よろしくお願い致します．

【main.m】

close all;

clc;

clear;

L = 10;

lim = L/2 - 0.5;

X = (combvec(0.1:0.1:0.5, -0.1:-0.1:-0.5, 0.5:0.5:lim, 0.5:0.5:lim, 0.5:0.5:lim, 0.5:0.5:lim, 0.5:0.5:lim))';

k1 = X(:, 1);

k2 = X(:, 2);

t1 = X(:, 3);

t2 = X(:, 4);

t3 = X(:, 5);

p1 = X(:, 6);

p2 = X(:, 7);

syms s k3

f2 = trap2(s, k1, k2, k3, t1, t2, t3, p1, p2, L);%%%%%% k3を求めるために変数をs,k3にした関数f2

eqn = int(f2, s, [0 L/2]) == 0; %%%%%%% k3に関する方程式を定義

k3 = solve(eqn,k3) ; %%%%%%%%%%%%%%%% 関数f2のトータル面積が0になるようなk3を求値

sprintf('fin')

【trap2.m】

function f = trap2(s, k1, k2, k3, t1, t2, t3, p1, p2, L)

t1_prime = t1 + p1;

t2_prime = t2 + p2;

a1 = (k2-k1)/(t2-t1_prime); %%% 傾き１

a2 = (k3-k2)/(t3-t2_prime); %%% 傾き２

f = piecewise(0<=s<=t1, (k1/t1)*s,...

t1<=s<=t1_prime, k1,...

t1_prime<=s<=t2,a1*(s-t1_prime)+k1,...

t2<=s<=t2_prime, k2,...

t2_prime<=s<=t3, a2*(s-t3)+k3,...

t3<=s<=L/2, k3);;

end

Answer 1

[Shota Kato]

以下の2点を確認させてください．

1. 台形の面積は解析的な解が存在するので，piecewise関数やint関数を用いることなく陽に答えを書くことができそうですが，そのようにしない回答が必要，ということでしょうか？そうでないなら，問題は簡単になります．
2. 添付画像を見ると，t1<t2<t3<L/2という条件が成立しそうですが，main.m内ではその条件が成り立っていないと思います．画像が正しいのであれば，コードを見直す必要がありますが，どうでしょうか？

[Shota Kato]

添付していただいたコードを回すと，確かにメモリ不足のエラーが出ますが，これはtrap2に与えられる引数がベクトルになっていて，関数内（15行目）の計算結果が行列になっているからです．

trap2の関数内の（/と*）の演算子をドット付き演算子にすることで，要素ごとに計算を行う（クロネッカー積と呼ばれます）ことができるので，メモリ不足の問題は解決されるはずです．

ただ，コードがこのままだとシンボリック計算が終了しない（主にintの部分）ため，計算がいつ終わるのかはわかりません．

[maro_ta]

ご回答ありがとうございます．ご質問を頂いたので回答させていただきます．

１，この先でtrap2の関数f(s)を用いた二重積分等の出力を行うのですが，この積分解としてパラメータを用いた解析解を与えるのは厳しいかと考え，台形を関数定義しております．

具体的には以下のような計算をするつもりでいます．

syms u

f = trap2(s, k1, k2, k3, t1, t2, t3, p1, p2, L); %%% k3求値後に関数fのを定義

theta = int(f, s, [0 u]);

x = double(int(cos(theta), u,[0 L]));

y = double(int(sin(theta), u,[0 L]));

２，おっしゃる通りですね．後々パラメータの大小関係を満たすものを抽出するつもりでいましたが，先に行った方が圧倒的にベクトルが短縮できたので，以下のように条件を満たすベクトルMを作りました(あまり得意でないので冗長になってしまいましたが....)

L = 10;

lim = L - 0.5;

X = (combvec(0.1:0.1:0.5, -0.1:-0.1:-0.5, 0.5:0.5:lim, 0.5:0.5:lim, 0.5:0.5:lim, 0.5:0.5:lim, 0.5:0.5:lim))'; %%% 順にk1,k2,t1,p1,t2,p2,t3

[Xrow Xcol] = size(X);

T = zeros(size(X));

count = 0;

for i=1:Xrow

if ((X(i,3) + X(i,4) < X(i,5))) && ((X(i,5) + X(i,6) < X(i,7))) %%% t1 + p1 < t2, t2 + p2 < t3 → ベクトルTを生成

count = count + 1;

T(count,:) = X(i,:);

end

end

clear X; %%% Xは不要なので削除

M = T(T(:,1) ~= 0,:); %%%% ベクトルTから要素が0の行を削除 → ベクトルMを生成

k1 = M(:, 1);

k2 = M(:, 2);

t1 = M(:, 3);

p1 = M(:, 4);

t2 = M(:, 5);

p2 = M(:, 6);

t3 = M(:, 7);

[maro_ta]

以下のコードで実行したところ，写真のようなエラーが発生しました．積分計算にベクトルを用いるのは厳しいのでしょうか....

【trapmat】

function f = trapmat(s, k1, k2, k3, t1, t2, t3, p1, p2, L)

t1_prime = t1 + p1;

t2_prime = t2 + p2;

a1 = (k2 - k1) ./ (t2 - t1_prime); %%% 傾き１

a2 = (k3 - k2) ./ (t3 - t2_prime); %%% 傾き２

f = piecewise(0<=s<=t1, (k1 ./ t1) .* s,...

t1<=s<=t1_prime, k1,...

t1_prime<=s<=t2, a1 .* (s - t2) + k2,...

t2<=s<=t2_prime, k2,...

t2_prime<=s<=t3, a2 .* (s - t3) + k3,...

t3<=s<=L/2, k3);

end

close all;
clc;
clear;

L = 10;
lim = L - 0.5;

X = (combvec(0.1:0.1:0.5, -0.1:-0.1:-0.5, 0.5:0.5:lim, 0.5:0.5:lim, 0.5:0.5:lim, 0.5:0.5:lim, 0.5:0.5:lim))'; %%% 順にk1,k2,t1,p1,t2,p2,t3
[Xrow Xcol] = size(X);
T = zeros(size(X)); 
count = 0;

for i=1:Xrow
    if ((X(i,3) + X(i,4) < X(i,5))) && ((X(i,5) + X(i,6) < X(i,7))) %%% t1 + p1 < t2, t2 + p2 < t3の条件でベクトルTを生成
        count = count + 1;
        T(count,:) = X(i,:);
    end
end

clear X; %%% Xは不要なので削除
M = T(T(:,1) ~= 0,:); %%%% ベクトルTから要素が0の行を削除 → ベクトルMを生成

k1 = M(:, 1);
k2 = M(:, 2);
t1 = M(:, 3);
p1 = M(:, 4);
t2 = M(:, 5);
p2 = M(:, 6);
t3 = M(:, 7);

syms s k3
f2 =  trapmat(s, k1, k2, k3, t1, t2, t3, p1, p2, L);%%%%%% k3を求めるために変数をs,k3にした関数f2
eqn = int(f2, s, [0 L/2]) == 0;  %%%%%%% k3に関する方程式を定義
k3 = solve(eqn,k3) ;  %%%%%%%%%%%%%%%% 関数f2のトータル面積が0になるようなk3を求値