4

私はmatlabを初めて使用し、近似x = aを開始してニュートンラフソン法をn回反復する関数を作成する必要があります。この開始近似は反復としてカウントされず、別の要件として for ループが必要です。投稿された他の同様の質問を見てきましたが、私の場合は while ループを使用したくありません。

これは私の入力が想定されているものです：

mynewton(f,a,n) which takes three inputs: 
f: A function handle for a function of x.
a: A real number.
n: A positive integer.

これまでのコードは次のとおりです。

function r=mynewton(f,a,n)
syms x;
z=f(x);
y=a;
for i=1:n    
    y(i+1)=y(i)-(z(i)/diff(z(i)));
end
r=y
end

関数を呼び出そうとすると、次のエラーが表示されます。

Error in MuPAD command: DOUBLE cannot convert the input expression into a double    array.
If the input expression contains a symbolic variable, use the VPA function instead.
Error in mynewton (line 6)
y(i+1)=y(i)-(z(i)/diff(z(i)));

問題は、この VPA 機能をどのように使用するかです。私のコードの残りの部分もおそらく100％正しいとは限りませんが、vpaの問題に対処したり、コードの他の部分を修正したりする助けがあれば大歓迎です.

ありがとう！

4

2 に答える 2

8

Newton-Raphson 法では正しくない点が 2 つありますが、確実に修正できます。これを修正した後VPA、あなたが話しているエラーは必要ありません。

エラー #1 - 反復更新

最初のものは反復そのものです。Newton-Raphson 法の定義を思い出してください。

何とか
(ソース: mit.edu )

次の反復では、前の反復の値を使用します。あなたがしていることは、ループ カウンターを使用してこれを に代入していることですがf(x)、これは正しくありません。これは、前の反復の値でなければなりません。

エラー #2 - 記号値と数値の混合

関数のコーディング方法を見ると、関数をシンボリックに定義しているにもかかわらず、数値を関数に代入しようとしています。残念ながら、これは MATLAB ではうまくいきません。実際に値を代入したい場合は、 を使用する必要がありますsubs。xこれにより、関数として、または関数が使用している独立変数の実際の値が置き換えられます。これを行うと、値はまだsymタイプです。これを数値的に使用できるようにするには、これをdoubleとしてキャストする必要があります。

yまた、効率のために、配列を作成する必要はありません。これを、反復ごとに更新される単一の値にするだけです。以上のことから、コードは次のように更新されます。念のために言っておきますが、必要な計算量を減らすために、ループの前に関数の導関数を取りました。また、ニュートン-ラフソン反復の分子項と分母項を分割して、物事を明確にし、これをsubs. 難しい話は抜きにして：

function r = mynewton(f,a,n)
syms x;
z = f(x);
diffZ = diff(z); %// Edit - Include derivative
y = a; %// Initial root

for idx = 1 : n    
    numZ = subs(z,x,y); %// Numerator - Substitute f(x) for f(y)
    denZ = subs(diffZ,x,y); %// Denominator - Substitute for f'(x) for f'(y)
    y = y - double(numZ)/double(denZ); %// Update - Cast to double to get the numerical value
end
r = y; %// Send to output
end

ループ内で置き換えiたことに注意してください。idxその理由は、これらの文字が複素数を表すために予約されているため、ループ インデックスとしてiorを使用することは実際には推奨されないためです。Shaijによるこの投稿を見ると、これらの変数をループ インデックスとして使用する方が実際には遅いことがわかります: Using i and j as variables in Matlab

いずれにせよ、これをテストするために、関数がy = sin(x)で、最初のルートが であると仮定しx0 = 2、5 回の反復で次のようにします。

f = @(x) sin(x);
r = mynewton(f, 2, 5)

r =

3.1416

sin(x)の切片は のsin(x)整数倍の位置にあるため、これは の知識と一致しますpi。 x0 = 2は近くにあるpiため、これは期待どおりに機能します。

あなたのための小さなボーナス

元のコードには、 の各反復でルートの値が格納されていましたy。本当にやりたい場合は、コードを次のように変更する必要があります。y物事をより効率的にするために事前に割り当てたことに注意してください。

function r = mynewton(f,a,n)
syms x;
z = f(x);
diffZ = diff(z);
y = zeros(1,n+1); %// Pre-allocate output array 
y(1) = a; %// First entry is the initial root

for idx = 1 : n    
    numZ = subs(z,x,y(idx)); %// Remember to use PREVIOUS guess for next guess
    denZ = subs(diffZ,x,y(idx));
    y(idx+1) = y(idx) - double(numZ)/double(denZ); %// Place next guess in right spot  
end
r = y; %// Send to output
end

上記とまったく同じパラメーターを使用してこのコードを実行すると、次のようになります。

f = @(x) sin(x);
r = mynewton(f, 2, 5)

r =

    2.0000    4.1850    2.4679    3.2662    3.1409    3.1416

の各値はr、その特定の反復におけるルートの推測を示します。配列の最初の要素は、最初の推測です (もちろん)。次の値は、ニュートン ラフソン ルートの各反復での推測です。配列の最後の要素が最後の反復であることに注意してください。これはほぼ に等しいpiです。

于 2014-08-05T03:23:15.263 に答える
0

使用する必要がありますeval

したがってz(x)=f(x)、eval(z(1))x = 1 で関数を評価するには

完全なコード:

function r=mynewton(fun,a,n)
% e.g. for mynewton(@sin, 1, 2)
syms x;
z(x)=fun(x);  % sin(x)
y=a;
for i=1:n    
    y(i+1)=y(i)-(eval(z(i))/eval(diff(z(i))));
end
r=y
end

質問の VPA 部分についてはよくわかりませんが、通常は MUPAD エラーを無視します:P

于 2014-08-05T01:50:25.250 に答える