matlab - quad を使用して関数を 1 つの変数だけで積分する

Question

2つの変数の関数を統合する方法はありますか?

f=@(x) x^2 + x*y

ちょうどxの上

quad(f, a, b) を試してみた

しかし、うまくいかず、代替ソリューションを探しています

Answer 1

次のようなものが必要なようです。

y = 100; % whatever y is
a = 0;
b = 2;
% you'll need to vectorize the integrand function
f = @(x) x.*x + x.*y
val = quad(f, a, b);

ただし、代数的な答えを探している場合は、Symbolic Toolbox、その他のソフトウェア、または微積分の本を使用する必要があります。:-)

「ベクトル化」すること全体は、Mathworksクワッドのドキュメントに基づいています。

関数y=fun（x）は、ベクトル引数xを受け入れ、ベクトル結果yを返す必要があります。これは、xの各要素で評価される被積分関数です。

Answer 2

申し訳ありませんが、quad は記号問題を解決しません。数値積分のみを行います。

syms x y
int(x^2 + x*y,x)
ans =
    (x^2*(2*x + 3*y))/6

シンボリックな問題を解決する自然な方法は、シンボリック ツールを使用することです。

フォローアップから、アーニャはその中間を望んでいます。ミックという古いロックスターの言葉を盗むために、「あなたはいつもあなたが望むものを手に入れることができるわけではありません.」

繰り返しますが、quad は適応ツールであるため、x のみで積分したい場合は quad を使用できません。

いくつかの単純なケースでは、シンプソンの法則のような単純なツールを使用して作業を行うことができます。たとえば、区間 [0 1] で x を積分して、上記の問題を解きたいとします。比較のために、まず象徴的に説明します。

syms x y
res = int(x^2 + x*y,x);

subs(res,x,1) - subs(res,0)
ans =
   y/2 + 1/3

では、x の数値積分を使って試してみましょう。

syms y
x = 0:.01:1;
coef = mod((0:100)',2)*2 + 2;
coef([1 end]) = 1;
coef = 0.01*coef/3;

(x.^2 + x.*y)*coef
ans =
   y/2 + 1/3

したがって、この非常に単純なケースでは、うまくいきました。もう少し複雑なものはどうですか？区間 [-1 1] で x*exp(x*y) を積分します。繰り返しますが、既知のフォームはシンボリックにアクセスできます。

syms x y
res = int(x*exp(x*y),x);

res = subs(res,x,1) - subs(res,-1)
res =
   (exp(-y)*(y + 1))/y^2 + (exp(y)*(y - 1))/y^2

後でテストするために、これは y = 1/2 でどのような値になりますか?

vpa(subs(res,y,1/2))
ans =
   0.34174141687554424792549563431876

シンプソンの法則を使って、同じトリックを試してみましょう。

syms y
x = -1:.01:1;
coef = mod((-100:100)',2)*2 + 2;
coef([1 end]) = 1;
coef = 0.01*coef/3;

res = (x.*exp(x*y))*coef
res =
   exp(y/2)/300 - exp(-y/2)/300 - exp(-y)/300 - exp(-y/4)/300 + exp(y/4)/300 - exp(-y/5)/750 + exp(y/5)/750 - exp(-(3*y)/4)/100 - exp(-(2*y)/5)/375 + exp((2*y)/5)/375 + exp((3*y)/4)/100 - exp(-(3*y)/5)/250 + exp((3*y)/5)/250 - (2*exp(-(4*y)/5))/375 + (2*exp((4*y)/5))/375 - exp(-y/10)/1500 + exp(y/10)/1500 - exp(-(3*y)/10)/500 + exp((3*y)/10)/500 - (7*exp(-(7*y)/10))/1500 + (7*exp((7*y)/10))/1500 - (3*exp(-(9*y)/10))/500 + (3*exp((9*y)/10))/500 - exp(-y/20)/1500 + exp(y/20)/1500 - exp(-(3*y)/20)/500 + exp((3*y)/20)/500 - exp(-y/25)/3750 + exp(y/25)/3750 - (7*exp(-(7*y)/20))/1500 - exp(-(2*y)/25)/1875 + exp((2*y)/25)/1875 + (7*exp((7*y)/20))/1500 - exp(-(3*y)/25)/1250 + exp((3*y)/25)/1250 - (3*exp(-(9*y)/20))/500 - (2*exp(-(4*y)/25))/1875 + (2*exp((4*y)/25))/1875 + (3*exp((9*y)/20))/500 - (11*exp(-(11*y)/20))/1500 - exp(-(6*y)/25)/625 + exp((6*y)/25)/625 + (11*exp((11*y)/20))/1500 - (7*exp(-(7*y)/25))/3750 + (7*exp((7*y)/25))/3750 - (13*exp(-(13*y)/20))/1500 - (4*exp(-(8*y)/25))/1875 + (4*exp((8*y)/25))/1875 + (13*exp((13*y)/20))/1500 - (3*exp(-(9*y)/25))/1250 + (3*exp((9*y)/25))/1250 - (11*exp(-(11*y)/25))/3750 + (11*exp((11*y)/25))/3750 - (17*exp(-(17*y)/20))/1500 - (2*exp(-(12*y)/25))/625 + (2*exp((12*y)/25))/625 + (17*exp((17*y)/20))/1500 - (13*exp(-(13*y)/25))/3750 + (13*exp((13*y)/25))/3750 - (19*exp(-(19*y)/20))/1500 - (7*exp(-(14*y)/25))/1875 + (7*exp((14*y)/25))/1875 + (19*exp((19*y)/20))/1500 - (8*exp(-(16*y)/25))/1875 + (8*exp((16*y)/25))/1875 - (17*exp(-(17*y)/25))/3750 + (17*exp((17*y)/25))/3750 - (3*exp(-(18*y)/25))/625 + (3*exp((18*y)/25))/625 - (19*exp(-(19*y)/25))/3750 + (19*exp((19*y)/25))/3750 - (7*exp(-(21*y)/25))/1250 + (7*exp((21*y)/25))/1250 - (11*exp(-(22*y)/25))/1875 + (11*exp((22*y)/25))/1875 - (23*exp(-(23*y)/25))/3750 + (23*exp((23*y)/25))/3750 - (4*exp(-(24*y)/25))/625 + (4*exp((24*y)/25))/625 - exp(-y/50)/7500 + exp(y/50)/7500 - exp(-(3*y)/50)/2500 + exp((3*y)/50)/2500 - (7*exp(-(7*y)/50))/7500 + (7*exp((7*y)/50))/7500 - (3*exp(-(9*y)/50))/2500 + (3*exp((9*y)/50))/2500 - (11*exp(-(11*y)/50))/7500 + (11*exp((11*y)/50))/7500 - (13*exp(-(13*y)/50))/7500 + (13*exp((13*y)/50))/7500 - (17*exp(-(17*y)/50))/7500 + (17*exp((17*y)/50))/7500 - (19*exp(-(19*y)/50))/7500 + (19*exp((19*y)/50))/7500 - (7*exp(-(21*y)/50))/2500 + (7*exp((21*y)/50))/2500 - (23*exp(-(23*y)/50))/7500 + (23*exp((23*y)/50))/7500 - (9*exp(-(27*y)/50))/2500 + (9*exp((27*y)/50))/2500 - (29*exp(-(29*y)/50))/7500 + (29*exp((29*y)/50))/7500 - (31*exp(-(31*y)/50))/7500 + (31*exp((31*y)/50))/7500 - (11*exp(-(33*y)/50))/2500 + (11*exp((33*y)/50))/2500 - (37*exp(-(37*y)/50))/7500 + (37*exp((37*y)/50))/7500 - (13*exp(-(39*y)/50))/2500 + (13*exp((39*y)/50))/2500 - (41*exp(-(41*y)/50))/7500 + (41*exp((41*y)/50))/7500 - (43*exp(-(43*y)/50))/7500 + (43*exp((43*y)/50))/7500 - (47*exp(-(47*y)/50))/7500 + (47*exp((47*y)/50))/7500 - (49*exp(-(49*y)/50))/7500 + (49*exp((49*y)/50))/7500 - exp(-y/100)/7500 + exp(y/100)/7500 - exp(-(3*y)/100)/2500 + exp((3*y)/100)/2500 - (7*exp(-(7*y)/100))/7500 + (7*exp((7*y)/100))/7500 - (3*exp(-(9*y)/100))/2500 + (3*exp((9*y)/100))/2500 - (11*exp(-(11*y)/100))/7500 + (11*exp((11*y)/100))/7500 - (13*exp(-(13*y)/100))/7500 + (13*exp((13*y)/100))/7500 - (17*exp(-(17*y)/100))/7500 + (17*exp((17*y)/100))/7500 - (19*exp(-(19*y)/100))/7500 + (19*exp((19*y)/100))/7500 - (7*exp(-(21*y)/100))/2500 + (7*exp((21*y)/100))/2500 - (23*exp(-(23*y)/100))/7500 + (23*exp((23*y)/100))/7500 - (9*exp(-(27*y)/100))/2500 + (9*exp((27*y)/100))/2500 - (29*exp(-(29*y)/100))/7500 + (29*exp((29*y)/100))/7500 - (31*exp(-(31*y)/100))/7500 + (31*exp((31*y)/100))/7500 - (11*exp(-(33*y)/100))/2500 + (11*exp((33*y)/100))/2500 - (37*exp(-(37*y)/100))/7500 + (37*exp((37*y)/100))/7500 - (13*exp(-(39*y)/100))/2500 + (13*exp((39*y)/100))/2500 - (41*exp(-(41*y)/100))/7500 + (41*exp((41*y)/100))/7500 - (43*exp(-(43*y)/100))/7500 + (43*exp((43*y)/100))/7500 - (47*exp(-(47*y)/100))/7500 + (47*exp((47*y)/100))/7500 - (49*exp(-(49*y)/100))/7500 + (49*exp((49*y)/100))/7500 - (17*exp(-(51*y)/100))/2500 + (17*exp((51*y)/100))/2500 - (53*exp(-(53*y)/100))/7500 + (53*exp((53*y)/100))/7500 - (19*exp(-(57*y)/100))/2500 + (19*exp((57*y)/100))/2500 - (59*exp(-(59*y)/100))/7500 + (59*exp((59*y)/100))/7500 - (61*exp(-(61*y)/100))/7500 + (61*exp((61*y)/100))/7500 - (21*exp(-(63*y)/100))/2500 + (21*exp((63*y)/100))/2500 - (67*exp(-(67*y)/100))/7500 + (67*exp((67*y)/100))/7500 - (23*exp(-(69*y)/100))/2500 + (23*exp((69*y)/100))/2500 - (71*exp(-(71*y)/100))/7500 + (71*exp((71*y)/100))/7500 - (73*exp(-(73*y)/100))/7500 + (73*exp((73*y)/100))/7500 - (77*exp(-(77*y)/100))/7500 + (77*exp((77*y)/100))/7500 - (79*exp(-(79*y)/100))/7500 + (79*exp((79*y)/100))/7500 - (27*exp(-(81*y)/100))/2500 + (27*exp((81*y)/100))/2500 - (83*exp(-(83*y)/100))/7500 + (83*exp((83*y)/100))/7500 - (29*exp(-(87*y)/100))/2500 + (29*exp((87*y)/100))/2500 - (89*exp(-(89*y)/100))/7500 + (89*exp((89*y)/100))/7500 - (91*exp(-(91*y)/100))/7500 + (91*exp((91*y)/100))/7500 - (31*exp(-(93*y)/100))/2500 + (31*exp((93*y)/100))/2500 - (97*exp(-(97*y)/100))/7500 + (97*exp((97*y)/100))/7500 - (33*exp(-(99*y)/100))/2500 + (33*exp((99*y)/100))/2500 + exp(y)/300

結果が得られましたが、それは私が望んでいた分析結果ではなく、少し厄介な混乱でした。それが正しいか？

vpa(subs(res,y,1/2))
ans =
   0.34174141693463006644516447861307

比較できるように、上記の分析結果をコピーします...

   0.34174141687554424792549563431876

ご覧のとおり、Simpson の規則は、[-1,1] に対して 0.01 のステップ サイズで、かなりうまく機能し、約 9 桁の 10 進数に一致しました。

この手法がより一般的なカーネルでも同様に機能するという保証はありませんが、必要なものが得られる可能性があります。