Finite difference method - solving boundary conditions

Hi,

I have problem writing finite differece method with next system of equations :

Does anyone have some ideas?

Thanks!

Answer 1

[Torsten]

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

[Maja Marevic]

Hi Torsten,

This works. But I want 3 graph (plots) for every function. Is it even possible with fsolve?

[Torsten]

plot(x,sol(1:n),x,sol(n+1:2*n),x,sol(2*n+1:3*n))