Hi. I'm attempting to solve the the ODEs one arrives at when analysing the problem of the flow due to a rotating disk, that is:
F''=F^2+G^2+F'H G''=2FG+G'H H''=-2F'
F(0)=0, G(0)=1, H(0)=0 F(infinity)=G(infinity)=0, H(infinity)=-A (constant)
I have implemented the BVP solver as such:
function M = Script(~)
M = bvpinit(linspace(0,4.5,801),@VKinit); sol = bvp4c(@VK,@VKbc,M);
function yprime = VK(~,y)
yprime = [ y(2) y(1)^2 - y(3)^2 + y(2)*y(5) y(4) 2*y(1)*y(3) + y(4)*y(5) y(6) -2*y(2)];
% y(1) = F, y(2) = F', y(3)= G, y(4) = G', y(5) = H and y(6)= H'.
function res = VKbc(ya,yb)
res = [ya(1);ya(3)-1;ya(5);yb(1);yb(3);yb(5)+0.88447];
function yinit = VKinit(~)
yinit = [0;0;0;0;0;0];
Plotting the solutions to this gives me exactly what I want, but that's only because I have taken the correct value of H(infinity)=-0.88447 from previous texts. (In this case I have taken infinity=4.5, the solution converges by this point)
I need to be able to find the value of this unknown boundary condition myself.
How do I go about modifying my code to obtain this value of H(infinty) for myself using a shooting method? Is the shooting method in fact the best way of evaluating H(infinity)?
Any help and advice anyone could give would be greatly appreciated!