finite difference method for second order ode

Hi everyone. I have written this code to solve this equation: y"+2y'+y=x^2 the problem is when I put X as for example X=0:0.25:1, it gives me fairly good answers for y. but when I change X as X=0:0.1:1, the answers for y are not correct. the more I reduce the delta x, the bigger the error become. can anyone tell me what I am doing wrong? this is the code:

%y"+2y'+y=x^2;

%Boundary Conditions:y(0)=0.2, y(1)=0.8;

x=0:0.25:1 %if I change 0.25 to 0.1 the answers are not acceptable.

n=length(x);

y=zeros(1,n);

y(1,1)=0.2; y(1,n)=0.8;

y

%((y(i+1)-2y(i)+y(i-1))/h^2)+2*(y(i+1)-y(i-1)/2h+y(i)=x(i)^2;

%after simplifing: 20y(i+1)-31y(i)+12y(i-1)=x(i)^2;

%AY=B ---->A=coefficients of y, Y=y's, B=the other side of the equations.

A=zeros(n-2);

B=zeros(1,n-2);

for i=1:n-2

    A(i,i)=-31;

end

for i=2:n-2

    A(i,i-1)=12;

    A(i-1,i)=20;

end

A %coefficient matrix

B(1,1)=x(1,2).^2-(12*y(1,1));

B(1,n-2)=x(1,n-1).^2-(20*y(1,n));

for i=2:n-3

    B(1,i)=(x(1,i+1)).^2;

end

B;

BB=B' %second side of the equations;

%AX=B ---> X=A\B

X=A\BB;

XX=X';

y(1,2:n-1)=XX(1,1:n-2);

y %final answers

plot(x,y,'-*');

Thank you in advance.

[madhan ravi]

Why do you say it’s incorrect?

[taher taher]

Thank you for the response. If you look at the pictures that I have attached, you can see the difference between the answers. I have solved the equation using "bvp4c" too and I know the answers should be like the first picture (h=0.25) also, does't reducing the delta x (h) mean that the answers should more precise?

[John D'Errico]

@Furkan Celik

Please do not answer questions just with a question. Use comments instead. Moved to a comment:

"Hi Taher!

Did you find the problem about it ?

I am trying to solve similar problem with your codes and ı am not sure is it correct or not.

Would you share script if it is possible...

Thank you..."

[Mustafa Ahmed]

Write MATLAB code to solve the following BVP using forward finite difference method:

𝑢′′ +1/𝑡 𝑢′ -1/𝑡^2 𝑢 = 0 𝑢(2) = 0.008

𝑢(6.5) = 0.003

ℎ = 1.5

And plot the approximate solution and the exact solution.

help here please

Answer 1

[John D'Errico]

Hmm. A zombie question, far too late for my answer to help the OP. But at least maybe I can help the secondary poster @Furkan Celik who revived the question from the long since dead. Odds are this will not really help though, since I cannot guess why Furcan had a problem. It might be unlikely exactly the same mistake was made twice.

This is a probably a case of sloppy programming, in the sense that Taher did not visualze the error being made. Of course, I don't have all of the code used. It is also a good reason not to write code to replace tools like an ODE solver, unless you have the programming skills and the numerical methods skills to have written that code yourself. And if you have those skills, you would also realize there is no reason to re-write that code in the first place.

But, lets look at what happened, where did the problem arise? I'd suggest it all comes from these two lines:

%((y(i+1)-2y(i)+y(i-1))/h^2)+2*(y(i+1)-y(i-1)/2h+y(i)=x(i)^2;
%after simplifing: 20y(i+1)-31y(i)+12y(i-1)=x(i)^2;

That is, the assumption that the discretized ODE does not change when you change the stride (h) in x. So the first line there had h in it. But the second line of code does not. Now, if you change the delta in x, things change.

Since this is a MATLAB forum and I am far too lazy to do actual arithmetic at this time of day, I'll use MATLAB to show the difference.

First, check the formula claimed by Taher:

syms y yp1 ym1
>> h = 0.25;

>> (yp1-2*y+ym1)/h^2 + 2*(yp1-ym1)/(2*h) + y
ans =
12*ym1 - 31*y + 20*yp1

So indeed the discretization redues to what was claimed. But, now if we change the stride, so a new h?

h = 0.1;

>> (yp1-2*y+ym1)/h^2 + 2*(yp1-ym1)/(2*h) + y
ans =
90*ym1 - 199*y + 110*yp1

So the coefficients in the discretization of the ODE are now different. If you used more elements in the vector x, but the OLD coefficients, you are essentially solving the wrong ODE.

Instead, better, more careful programming practice would not have allowed this mistake. By hard coding those coefficients into the code as we see was done by Taher:

for i=1:n-2
    A(i,i)=-31;
end
for i=2:n-2
    A(i,i-1)=12;
    A(i-1,i)=20;
end
A %coefficient matrix
B(1,1)=x(1,2).^2-(12*y(1,1));
B(1,n-2)=x(1,n-1).^2-(20*y(1,n));

you are virtually asking for troubles, asking to create buggy code. I'm sorry, but that is just code pleading for bugs to arise.

So, now let me quickly slap together a code to solve this problem in a way that will easily allow me to change the stride. Then we can see if it works for various step sizes. I'll bet it does.

I'll first write a function to solve this ODE using an indicated discretization, since the solution will need to be done several times, for various step sizes.

function [y,x] = mybvpsolve(xlim,y0,y1,n)
  % y"+2y'+y=x^2
  % xlim is a vector of length 2
  % y0, y1 are the BC at each end
  % n is the number of elements in the final result
  
  x = linspace(xlim(1),xlim(2),n)';
  
  % stride
  h = x(2) - x(1);
  
  % discretize the ODE into a tridiagonal matrix
  % ((y(i+1)-2y(i)+y(i-1))/h^2) + 2*(y(i+1)-y(i-1)/2h + y(i)=x(i)^2;
  
  % so the main diagonal has coefficients
  % the 1 at each end brings the boundary conditions into the problem.
  maindiag = [1;repmat(-2/h^2 + 1,n - 2,1);1];
  % be careful to create the sub and super diagonals as the correct lengths
  subdiag = [repmat(1/h^2 - 1/h,n-2,1);0];
  supdiag = [0;repmat(1/h^2 + 1/h,n-2,1)];
  
  % 3 calls to diag will do the trick here. For large problems, a
  % sparse matrix would be faster and better yet.
  ODEmat = diag(maindiag) + diag(subdiag,-1) + diag(supdiag,1);
  
  % create the right hand side for the solve:
  rhs = [y0;x(2:end-1).^2;y1];
  
  % solve. This is just a call to backslash
  y =ODEmat\rhs;
  
end

Now, call it 3 times, with various numbers of nodes in the solution, and then plot the results.

xlim = [0 1];
y0 = 0.2;
y1 = 0.8;

% 0:.25:1
[y5,x5] = mybvpsolve(xlim,y0,y1,5);
% 0:.1:1
[y11,x11] = mybvpsolve(xlim,y0,y1,11);
% 0:.01:1
[y101,x101] = mybvpsolve(xlim,y0,y1,101);

plot(x5,y5,'r-',x11,y11,'go',x101,y101,'b.')

What you see here is that the solution with 11 nodes is almost indistinguishable from the solution with 101 nodes. Even the red solution curve with only 5 nodes did quite well. And all are on top of each other.

In the end, I strongly recommend you just use the canned BVP4c solver.

[Furkan Celik]

Thank you so much for your response

[kibuchi]

Thank you for this

.