MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn moreOpportunities for recent engineering grads.

Apply Today**New to MATLAB?**

Asked by pxg882
on 22 Oct 2012

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'

with BCs

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.*

## 0 Comments