Finite difference method - solving boundary conditions

11 ビュー (過去 30 日間)
Maja Marevic
Maja Marevic 2019 年 7 月 8 日
コメント済み: Torsten 2019 年 7 月 9 日
Hi,
I have problem writing finite differece method with next system of equations :
Does anyone have some ideas?
Thanks!

採用された回答

Torsten
Torsten 2019 年 7 月 8 日
編集済み: Torsten 2019 年 7 月 8 日
y1(1) - 1.0 = 0
(y1(i+1)-2*y1(i)+y1(i-1))/dx^2 - (2*(y1(i+1)-y1(i-1))/(2*dx)/(1-x(i)) + y1(i)/(x(i)^2*(1-x(i))^2)*(y1(i)^2-1+(y3(i)^2-y2(i)^2)/(1-x(i))^2)) = 0
(i = 2,...,n-1)
y1(n) = 0
y2(1) = 0
(y2(i+1)-2*y2(i)+y2(i-1))/dx^2 - (2*y2(i)*y1(i)^2/(x(i)^2*(1-x(i))^2)) = 0
(i=2,...,n-1)
y2(n) - eta = 0
y3(1) = 0
(y3(i+1)-2*y3(i)+y3(i-1))/dx^2 - (y3(i)/(x(i)^2*(1-x(i))^2)*(2*y1(i)^2+beta*(y3(i)^2-x(i)^2)/(1-x(i))^2)) = 0
(i=2,...,n-1)
y3(n) - 1.0 = 0
with
x(i) = (i-1)/(n-1) (i=1,...,n)
3*n nonlinear equations in 3*n unknowns.
You can try "fsolve" to solve:
function main
n = 100;
x0 = ones(3*n,1);
sol = fsolve(@(x)fun(x,n),x0);
norm(fun(sol,n))
x = ((1:n)-1)/(n-1);
plot(x,sol(1:n))
end
function res = fun(z,n)
eta=1.0;
beta = 1.0;
x = ((1:n)-1)/(n-1);
dx = 1/(n-1);
y1 = z(1:n);
y2 = z(n+1:2*n);
y3 = z(2*n+1:3*n);
res_y1 = zeros(n,1);
res_y2 = zeros(n,1);
res_y3 = zeros(n,1);
res_y1(1) = y1(1)-1.0;
for i=2:n-1
res_y1(i) = (y1(i+1)-2*y1(i)+y1(i-1))/dx^2 - (2*(y1(i+1)-y1(i-1))/(2*dx)/(1-x(i)) + y1(i)/(x(i)^2*(1-x(i))^2)*(y1(i)^2-1+(y3(i)^2-y2(i)^2)/(1-x(i))^2));
end
res_y1(n) = y1(n);
res_y2(1) = y2(1);
for i = 2:n-1
res_y2(i) = (y2(i+1)-2*y2(i)+y2(i-1))/dx^2 - (2*y2(i)*y1(i)^2/(x(i)^2*(1-x(i))^2));
end
res_y2(n) = y2(n)-eta;
res_y3(1) = y3(1);
for i=2:n-1
res_y3(i) = (y3(i+1)-2*y3(i)+y3(i-1))/dx^2 - (y3(i)/(x(i)^2*(1-x(i))^2)*(2*y1(i)^2+beta*(y3(i)^2-x(i)^2)/(1-x(i))^2));
end
res_y3(n) = y3(n)-1.0;
res = [res_y1;res_y2;res_y3];
end
  2 件のコメント
Maja Marevic
Maja Marevic 2019 年 7 月 9 日
Hi Torsten,
This works. But I want 3 graph (plots) for every function. Is it even possible with fsolve?
Torsten
Torsten 2019 年 7 月 9 日
plot(x,sol(1:n),x,sol(n+1:2*n),x,sol(2*n+1:3*n))

サインインしてコメントする。

その他の回答 (0 件)

カテゴリ

Help Center および File ExchangeLinear Algebra についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by