MATLAB Answers

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

3 views (last 30 days)
maro_ta
maro_ta on 11 Oct 2020
Commented: maro_ta on 12 Oct 2020
Xにより7つのパラメータを遷移させた上で,解としてk3を求めるプログラムを組んだのですが,エラーが発生してしまいました.
操作としては,Xで全てのパラメータのパターンを網羅したベクトルをつくり,別のmファイル(trap2.m)から関数を呼び出し各パラメータを代入し,方程式eqnを満たすk3を求めるといったものになります.
また,方程式eqnはパラメータk1,k2,t1,t2,t3,p1,p2を与えた上で,trap2.m(台形を組み合わせた関数)の面積が0となるような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); %%% 傾き1
a2 = (k3-k2)/(t3-t2_prime); %%% 傾き2
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

Answers (1)

Shota Kato
Shota Kato on 11 Oct 2020
以下の2点を確認させてください.
  1. 台形の面積は解析的な解が存在するので,piecewise関数やint関数を用いることなく陽に答えを書くことができそうですが,そのようにしない回答が必要,ということでしょうか?そうでないなら,問題は簡単になります.
  2. 添付画像を見ると,t1<t2<t3<L/2という条件が成立しそうですが,main.m内ではその条件が成り立っていないと思います.画像が正しいのであれば,コードを見直す必要がありますが,どうでしょうか?
  3 Comments
maro_ta
maro_ta on 12 Oct 2020
以下のコードで実行したところ,写真のようなエラーが発生しました.積分計算にベクトルを用いるのは厳しいのでしょうか....
【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); %%% 傾き1
a2 = (k3 - k2) ./ (t3 - t2_prime); %%% 傾き2
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を求値

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!