1

私たちは、MATLAB で独自の微分演算子を定義するように求められました。私は一連の手順に従ってそれを行いました。次に、微分演算子を使用して境界値の問題を解決する必要があります。

-y'' + 2y' - y = x、y(0) = y(1) =0

私のコードは次のとおりで、これを計算するために使用されました（一次および二次導関数）

h = 2;
x = 2:h:50; 
y = x.^2 ;

n=length(x);
uppershift = 1; 
U = diag(ones(n-abs(uppershift),1),uppershift);
lowershift = -1;
L= diag(ones(n-abs(lowershift),1),lowershift);
% the code above creates the upper and lower shift matrix


D  = ((U-L))/(2*h); %first differential operator 
D2 = (full (gallery('tridiag',n)))/ -(h^2); %second differential operator

d1= D*y.'
d2= ((D2)*y.')

ここに投稿した後、これに変更し、Identity Matrix の使用を奨励する 1 つの応答を得ましたが、まだどこにも到達していないようです。

h = 2;

n=10;
uppershift = 1; 
U = diag(ones(n-abs(uppershift),1),uppershift);
lowershift = -1;
L= diag(ones(n-abs(lowershift),1),lowershift);


D  = ((U-L))/(2*h); %first differential operator 
D2 = (full (gallery('tridiag',n)))/ -(h^2); %second differential operator


I= eye(n);

eqn=(-D2 + 2*D - I)*y == x

solve(eqn,y)

y と x を定義する必要があるか、正確には何を定義する必要があるかなど、これをどのように進めればよいかわかりません。私は無知です！

4

1 に答える 1

1

これはODE の解の数値x=0近似であるため、時間からまでのこの ODE の解を表す数値ベクトルを見つけようとしていx=1ます。これは、境界条件によって、解が 0 と 1 の間でのみ有効になることを意味します。

また、これは現在逆の問題です。 一緒に行った前回の投稿では、入力ベクトルが何であるかを知っており、行列とベクトルの乗算を行うと、その入力ベクトルで出力導関数演算が生成されました。ここで、導関数の出力が与えられ、元の入力が何であったかを求めています。これには、線形連立方程式を解くことが含まれます。

基本的に、あなたの問題は次のとおりです。

YX = F

Yは、導出した行列導関数演算子からの係数です。これはn x n行列でXあり、ODE の解になります。これはn x 1ベクトルでありF、ODE を関連付ける関数でもあり、n x 1ベクトルでもあります。私たちの場合、それはx. を見つけるYには、コード内ですでにほとんどのことを行っています。各行列演算子 (一次導関数と二次導関数) を単純に取り、ODE の左辺を尊重する適切な符号とスケールと共にそれらを加算します。ところで、一次導関数と二次導関数の行列は正しいです。-y残っているのは、用語をミックスに追加すること-eye(n)です。これは、コードで見つけたように達成されます。

Yandを定式化したら、 or演算子をF使用してを解き、次の方法でこの線形システムの解を得ることができます。mldivide\X

X = Y \ F;

上記は基本的に、 と によって形成される方程式の線形システムを解き、YにF格納されXます。

x=0最初に行う必要があるのは、 からまでの点のベクトルを定義することですx=1。 linspace必要なポイント数を指定できる場合は、おそらく最も適しています。とりあえず100点としましょう。

x = linspace(0,1,100);

したがって、h私たちの場合は です1/100。x = a一般に、始点 から終点まで解きたい場合x = b、ステップ サイズhは次のように定義されますh = (b - a)/n。ここnで、 は ODE で解きたい点の総数です。

ここで、境界条件を含める必要があります。これは単純に、ODE の解の始まりと終わりを知っていることを意味します。これは、y(0) = y(1) = 0. そのため、 の最初の行のY最初の列のみが 1 に設定され、 の最後の行の最後の列Yのみが 1 に設定されていることを確認し、両方の出力位置をF0 に設定します。これは、これらの点での解決策はすでにわかっています。

したがって、解決する最終的なコードは次のとおりです。

%// Setup
a = 0; b = 1; n = 100;
x = linspace(a,b,n);
h = (b-a)/n;

%// Your code
uppershift = 1;
U = diag(ones(n-abs(uppershift),1),uppershift);
lowershift = -1;
L= diag(ones(n-abs(lowershift),1),lowershift);
D  = ((U-L))/(2*h); %first differential operator
D2 = (full (gallery('tridiag',n)))/ -(h^2);

%// New code - Create differential equation matrix
Y = (-D2 + 2*D - eye(n));

%// Set boundary conditions on system
Y(1,:) = 0; Y(1,1) = 1;
Y(end,:) = 0; Y(end,end) = 1;

%// New code - Create F vector and set boundary conditions
F = x.';
F(1) = 0; F(end) = 0;

%// Solve system
X = Y \ F;

Xh = 1/100からx=0までのステップで、ODE の数値近似が含まれているはずx=1です。

これがどのように見えるか見てみましょう：

figure; 
plot(x, X);
title('Solution to ODE');
xlabel('x'); ylabel('y');

ここに画像の説明を入力

y(0) = y(1) = 0境界条件に従ってそれを見ることができます。

これがお役に立てば幸いです。幸運を祈ります。

于 2015-06-28T19:29:33.793 に答える