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

1 次查看（过去 30 天）
maro_ta
maro_ta 2020-10-11
评论： maro_ta 2020-10-12
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

请先登录，再回答此问题。

回答（1 个）

Shota Kato
Shota Kato 2020-10-11
以下の2点を確認させてください．
  1. 台形の面積は解析的な解が存在するので，piecewise関数やint関数を用いることなく陽に答えを書くことができそうですが，そのようにしない回答が必要，ということでしょうか？そうでないなら，問題は簡単になります．
  2. 添付画像を見ると，t1<t2<t3<L/2という条件が成立しそうですが，main.m内ではその条件が成り立っていないと思います．画像が正しいのであれば，コードを見直す必要がありますが，どうでしょうか？
  3 个评论
显示 1更早的评论
maro_ta
maro_ta 2020-10-12
ご回答ありがとうございます．ご質問を頂いたので回答させていただきます．
１，この先で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
maro_ta 2020-10-12
以下のコードで実行したところ，写真のようなエラーが発生しました．積分計算にベクトルを用いるのは厳しいのでしょうか....
【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を求値

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 微積分 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!