Plotting graphs from equation of bvp

Hi, I have some problems regarding bvp. It's all about plotting graphs.

This is my code:

function Script

n=0.6;
infinity=20;

solinit = bvpinit(linspace(0,infinity,30),[0 0 0 0 -1]);

sol = bvp4c(@(x,y)VK(x,y,n),@VKbc,solinit);
x = [0:2:20];
y = deval(sol,x);
figure
title('Plots of U versus x') 
ylabel('U(x)')
xlabel('x')
drawnow
hold on

for i=1:9
  n = n+0.1;
  sol = bvp4c(@(x,y)VK(x,y,n),@VKbc,sol);
  plot(sol.x,sol.y(1,:));

    end

hold off

function yprime = VK(x,y,n)
a = ((1-n)/(n+1))*x;
c = (y(4)^2+y(5)^2)^((1-n)/(n+1));
yprime = [ c*y(4)
           c*y(5)
           -2*y(1) - a*c*y(4)
           y(1)^2 - (y(2)+1)^2 + (y(3)+a*y(1))*c*y(4)
           2*y(1)*(y(2)+1) + (y(3)+a*y(1))*c*y(5)];
end

function res = VKbc(ya,yb)

res = [ya(1);ya(2);ya(3);yb(1);yb(2)+1];

end

end

My problem is:

1) I have to plot c aginst x.

Where:

c = (y(4)^2+y(5)^2)^((1-n)/(n+1))

How should I write the code? 2) I also have to plot U', V' and W', all aginst x. Where, (as from the code)

y(1)'= U'=c*y(4),
y(2)' = V'=c*y(5),
y(3)' = W'=-2*y(1) - a*c*y(4).

I hope someone can help me regarding this. Thank you in advance.

Answer 1

[Star Strider]

This works:

function Script

n=0.6;
infinity=20;

solinit = bvpinit(linspace(0,infinity,30),[0 0 0 0 -1]);

sol = bvp4c(@(x,y)VK(x,y,n),@VKbc,solinit);
x = [0:2:20];
y = deval(sol,x);

for i=1:9
    n = n+0.1;
    sol = bvp4c(@(x,y)VK(x,y,n),@VKbc,sol);

      xc{i} = sol.x;
      yc{i} = sol.y;

      cv{i} = (sol.y(4,:).^2+sol.y(5,:).^2).^((1-n)./(n+1));          % Calculate ‘c’

      figure(i)
      plot(xc{i},yc{i}, xc{i},cv{i});
      title('Plots of U versus x')
      ylabel('U(x)')
      xlabel('x')
      legend('y_1(x)', 'y_2(x)', 'y_3(x)', 'y_4(x)', 'y_5(x)', 'c(x)', 'Location','NE')

end

      function yprime = VK(x,y,n)
      a = ((1-n)/(n+1))*x;
      c = (y(4)^2+y(5)^2)^((1-n)/(n+1));
      yprime = [ c*y(4)
          c*y(5)
          -2*y(1) - a*c*y(4)
          y(1)^2 - (y(2)+1)^2 + (y(3)+a*y(1))*c*y(4)
          2*y(1)*(y(2)+1) + (y(3)+a*y(1))*c*y(5)];
      end

      function res = VKbc(ya,yb)

      res = [ya(1);ya(2);ya(3);yb(1);yb(2)+1];

      end
  end

I calculated ‘c’, and created cell arrays to save the solutions at each value of ‘i’. You can use ‘xc’, ‘yc’, and ‘cv’ in other parts of your code. See the documentation on cell arrays for details on how to work with them. (In this instance, they should be relatively easy to work with. The plot call in the loop is an example of how to address them.)

I don’t know how you want your solutions to be plotted. Here, I plotted 9 separate figures, since it would be extremely difficult for you to see the individual plots if you plotted them all in one figure. I also added a legend call to each figure so you could see what the lines are.

[Day Rosli]

Thankyou so much! Actually, I want to plot for different number of n. Eg: Graph of U against x with n=0.6-1.4.

But, this really help me alot though. Thank you again!

[Star Strider]

My pleasure!

If you want to plot ‘U’ as a function of ‘n’ and ‘x’, it is likely best to use the ribbon plot.

[Day Rosli]

Oh, unfortunately this is my first time hearing about ribbon plot. As for the time being, since I am only used to the basic plot, I think I am gonna plot the U', V' and W' just like what you teach me on how to plot c. Since

U'=c*y(4),
V'=c*y(5),
W'=-2*y(1) - a*c*y(4).

So, it will be like;

plot(sol.x,(sol.y(4,:).^2+sol.y(5,:).^2).^((1-n)./(n+1)).*sol.y(4,:));

I know there must be other better way but since this is my first time using MATLAB, so I will just do like this for now.

Thank you for your idea and time! :)

[Star Strider]

My pleasure!

To calculate U',V',W', your loop changes to:

for i=1:9
    n = n+0.1;
    sol = bvp4c(@(x,y)VK(x,y,n),@VKbc,sol);

      xc{i} = sol.x;
      yc{i} = sol.y;
      L{i} = length(xc{i});

      cv{i} = (sol.y(4,:).^2+sol.y(5,:).^2).^((1-n)./(n+1));          % Calculate ‘c’
      Udot{i} = cv{i} .* yc{i}(4,:);
      Vdot{i} = cv{i} .* yc{i}(5,:);
      av = ((1-n)./(n+1))*xc{i};
      Wdot{i} = -2*yc{i}(1,:) - av.*cv{i}.*yc{i}(4,:);

      figure(i)
      plot(xc{i},yc{i}, xc{i},cv{i});
      title('Plots of U versus x')
      ylabel('U(x)')
      xlabel('x')
      legend('y_1(x)', 'y_2(x)', 'y_3(x)', 'y_4(x)', 'y_5(x)', 'c(x)', 'Location','NE')

end

Where ‘Udot’ = U', and so for the others. (The rest of your code is unchanged, and with these changes, you complete code ran for me without error.) I added the ‘L’ vector (the length of each ‘x’ vector), since you will need that information if you want to create a double array to do a ribbon or other 3D plot. (I did not separately save ‘av’, since you do not seem to need it other than to calculate ‘Wdot’.)

The statements appear correct to me, however you need to check them to be sure I wrote them correctly. You can plot them by their names, since the vectors (arrays for different values of ‘n’) now exist. If everything is calculated correctly, I would save all the variables you want to use to a .mat file so you do not have to re-run your bvp4c code each time. See the documentation on save, load, and matfile for all the interesting details.

I don’t know what you’re doing, but I’m sure you’ll succeed with it. It looks interesting.

[Day Rosli]

Again, thank you for your help! Your explanation really help me to write this program.

By the way, as I want to add one parameter which is M in the equation. And as I vary the values of M, I got error like

Unable to solve the collocation equations -- a singular Jacobian encountered.

Mind to explain to me what is actually going on?

In order to stop from getting this error, I tried to change the value of the parameter to other values. But still, if only I know how to deal with this kind of error, I could just use the values which I should use at the first place.

For example:

I need to use values of M (which I need to add into the equation above) 
for M=1,4,7,10. But i got the singular jacobian error. 
So instead, I choose M=0.1,0.2,0.3,0.4. And it works well without error. 
I can just proceed with the smaller values of M but if only I know why and 
how to correct my code without getting such error, it would be great!

Thank you :)

[Star Strider]

I’m not certain, but since it depends on the magnitude of ‘M’, I assume that the current Jacobian has a very high condition number, and adding ‘M’ with large enough values is enough for it to lose rank. I haven’t delved into the code for bvp4c, so I can’t say for sure. I would just go with the smaller values of ‘M’. (Computers have their limits.)

If you want to experiment with the Jacobian on your own in detail, and if you have the Symbolic Math Toolbox (that makes the calculations easier and does not create algebra errors), you can calculate the symbolic Jacobian, and using either the symbolic or numeric version (use matlabFunction to get it as a function you can use in core MATLAB), explore its properties with various values of your constants. Use the cond function to get the condition number. Lower condition numbers are better.

My pleasure.

[Day Rosli]

Thank you so much again, Star Strider. I believe, you have share alot of your knowledge to others. And Im one of them and Im so glad. Wish you the best in your life and career!

:)

[Star Strider]

As always, my pleasure!

Thank you! I wish you the same!