Hi below is matlab code for the boundary value problem y'''' + yy'' + x2 cos (pi(y)) -1 = 0; y(-1) = 0; y(1) = 0; y'(-1) = 1; y'(1) = 1:
I want the code to display a matrix SOL whose columns are successive approximations
to the solution. But the matrix just isn't showing up...can anyone see why.
Thanks
% This is a boilerplate code to calculate the solution to a
% nonlinear BVP of fourth order, two dirichlet and two neumann BCs
% Define the interval over which solution is calculated
a= -1;
b= 1;
% Define N
N=100;
% Define the grid spacing
h=(b-a)/N;
% Define the grid
x= reshape ( linspace (a,b,N+1) ,N+1 ,1);
% Define an initial guess for the solution of the BVP (make sure it
% is a column vector)
U=x.*(1-x);
% Store the initial guess in SOL
SOL =U;
% Define a tolerance for the termination of Newton?Raphson
tol =1e-4;
% Ensure that F is such that we always do at least one iteration
F= ones (N+1 ,1);J= zeros (N+1,N +1);
% The loop continues while ever the norm of F is greater than tol
while ( norm (F)>tol)
%% Define F
% Boundary conditions
F (1)= U(1) ;
F(N +1)=U(N+1);
% Boundary conditions (Fictitious nodes) + ODE
F (2)=((3*U(1)-4*U(2)-U(3)+2*h)/h^2) - (1-x(1));
F(N)=((3*U(N+1)-4*U(N)+U(N-1)-2*h)/h^2)-(1-x(1));
% Finite difference approximation to ODE at interior nodes
for i =3N -1),
F(i)=(U(i+2)-4*U(i+1)+6*U(i)-4*U(i-1)+U(i-2))+(h^2*U(i))*(U(i+1)-2*U(i)+U(i-1))+(h^4*x(i))*(cos(pi*U(i)))-(h^4);
end
%% Define J
% First row corresponds to F(1)=...;
J (1 ,1)=1;
% Last row corresponds to F(N+1)=...;
J(N +1 ,N+1)=1;
% Second row corresponds to F(2)=...;
J (2 ,1)=3/h^2;
J (2 ,2)=-4/h^2;
J (2 ,3)=-1/h^2;
% Penultimate row corresponds to F(2)=...;
J(N ,N-1)=1/h^2;
J(N ,N)=-4/h^2;
J(N ,N+1)=3/h^2;
% Intermediate rows correspond to F(i)=...;
for i =3N -1) ,
% Diagonal entry
J(i,i)=6-2*h^2*U(i)^2 +h^4*x(i)*cos(pi*U(i));
J(i,i-2)=1;
J(i,i-1)=-4+h^2*U(i)*U(i-1);
J(i,i+1)=-4+h^2*U(i)*U(i+1);
J(i,i+2)=1;
% There may be off?diagonal entries, check your calculation
end
% Having defined F and J, update the approximate solution to
% difference equations
U=U-J\F;
% Store the new approximation
SOL =[ SOL ,U];
end
plot(x,U);
% Insert your plot command
I want the code to display a matrix SOL whose columns are successive approximations
to the solution. But the matrix just isn't showing up...can anyone see why.
Thanks
% This is a boilerplate code to calculate the solution to a
% nonlinear BVP of fourth order, two dirichlet and two neumann BCs
% Define the interval over which solution is calculated
a= -1;
b= 1;
% Define N
N=100;
% Define the grid spacing
h=(b-a)/N;
% Define the grid
x= reshape ( linspace (a,b,N+1) ,N+1 ,1);
% Define an initial guess for the solution of the BVP (make sure it
% is a column vector)
U=x.*(1-x);
% Store the initial guess in SOL
SOL =U;
% Define a tolerance for the termination of Newton?Raphson
tol =1e-4;
% Ensure that F is such that we always do at least one iteration
F= ones (N+1 ,1);J= zeros (N+1,N +1);
% The loop continues while ever the norm of F is greater than tol
while ( norm (F)>tol)
%% Define F
% Boundary conditions
F (1)= U(1) ;
F(N +1)=U(N+1);
% Boundary conditions (Fictitious nodes) + ODE
F (2)=((3*U(1)-4*U(2)-U(3)+2*h)/h^2) - (1-x(1));
F(N)=((3*U(N+1)-4*U(N)+U(N-1)-2*h)/h^2)-(1-x(1));
% Finite difference approximation to ODE at interior nodes
for i =3N -1),
F(i)=(U(i+2)-4*U(i+1)+6*U(i)-4*U(i-1)+U(i-2))+(h^2*U(i))*(U(i+1)-2*U(i)+U(i-1))+(h^4*x(i))*(cos(pi*U(i)))-(h^4);
end
%% Define J
% First row corresponds to F(1)=...;
J (1 ,1)=1;
% Last row corresponds to F(N+1)=...;
J(N +1 ,N+1)=1;
% Second row corresponds to F(2)=...;
J (2 ,1)=3/h^2;
J (2 ,2)=-4/h^2;
J (2 ,3)=-1/h^2;
% Penultimate row corresponds to F(2)=...;
J(N ,N-1)=1/h^2;
J(N ,N)=-4/h^2;
J(N ,N+1)=3/h^2;
% Intermediate rows correspond to F(i)=...;
for i =3N -1) ,
% Diagonal entry
J(i,i)=6-2*h^2*U(i)^2 +h^4*x(i)*cos(pi*U(i));
J(i,i-2)=1;
J(i,i-1)=-4+h^2*U(i)*U(i-1);
J(i,i+1)=-4+h^2*U(i)*U(i+1);
J(i,i+2)=1;
% There may be off?diagonal entries, check your calculation
end
% Having defined F and J, update the approximate solution to
% difference equations
U=U-J\F;
% Store the new approximation
SOL =[ SOL ,U];
end
plot(x,U);
% Insert your plot command