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!
No products are associated with this question.