MATLAB: Finite difference method – solving boundary conditions

bvpdifferential equationsfinite difference methodsystem of differential equations

Hi,
I have problem writing finite differece method with next system of equations :
Does anyone have some ideas?
Thanks!

Best Answer

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
Related Question