Thread: Shooting method

hi guys,
i must di this exercise:
Solve with shooting method this boundary value problem:

[V-(1-m)*(y/2)]*(dU/dy)+m*U^2-d/dy(dU/dy)=0
-[(1-m)*(y/2)]*(dU/dy)+m*U+dV/dy=0

for U(y) and V(y) with y[0..20]

boundary value: U(0)=0, V(0)=0 and U(20)=1

with

1 - Euler method
2 - Rounge Kutta 4 order


function [u v w y]=shooting_Eulero_Esp(dy)
global u v w y L N h uL m


h=0.01;


L=20;


N=(L/h)+1;


y=zeros(N,1);
for i=1:length(y)


y(i)=h*(i-1);
end


u=zeros(N,1);


v=zeros(N,1);
w=zeros(N,1);
m=0;
toll=10^-12;


u(1)=0;
v(1)=0;
uL=1;


eta1=1;
eta2=2;


[check1 err1]=shoot_EE(eta1);


[check2 err2]=shoot_EE(eta2);

uN=u(L)


numiter=0;
err=1;


while(abs(err)>=toll)
numiter=numiter+1;
etat=eta2-err2*(eta1-eta2)/(err1-err2);


[check err]=shoot_EE(etat);
if(abs(err)<abs(err2))
eta1=eta2;
err1=err2;
eta2=etat;
err2=err;
else
eta1=etat;
err1=err;
end
end
numiter
uN=u(length(u))
w0=w(1)
plot(y,u,y,w)
axis([0 20 0 1.2])



function [uN err]=shoot_EE(eta)
global y u v w h uL m


w(1)=eta;


for i=2:length(u)
u(i)=u(i-1)+h*w(i-1);
v(i)=v(i-1)+h*(((1-m)*y(i-1)*w(i-1)/2)-m*u(i-1));
w(i)=w(i-1)+h*(v(i-1)*w(i-1)-((1-m)*u(i-1)*w(i-1)*y(i-1)/2)+m*(u(i-1)^2)+m);
end
uN=u(length(u));
err=uL-uN;
abs(err)
return


CAN SOMEONE HELP ME WITH RUNGE KUTTA 4 ORDER?

Hey nunzio.

Is this in MATLAB? If so did you plot the results of Euler and RK4? If Euler works then is the RK4 resembling the Euler plot?

hi chiro,
yes this is matlab code. in this exercise i need to compare the solution of U and dU/dy for both method.
for euler method i have this MATLAB.pdf

i need to implement a runge kutta 4 order

thanks

Can you point out the routine in pseudo-code? (For RK4 that is).

i am sorry but i dont have tue routine for rk4

I dug up a MATLAB routine from my old course-work:

k1 = h*f(x[j],y[j]);
k2 = h*f(x[j]+0.5*h,y[j] + 0.5*k1);
k3 = h*f(x[j] + 0.5*h,y[j] + 0.5*k2);
k4 = h*f(x[j] + h, y[j] + k3);
y[j+1] = y[j] + (1/6)*(k1 + 2*k2 + 2*k3 + k4);

In this example, h is the step-size x[j] the array of x values and y[j] the function values for corresponding x[j].

in my exercise there are
[V-(1-m)*(y/2)]*(dU/dy)+m*U^2-d/dy(dU/dy)=0
-[(1-m)*(y/2)]*(dU/dy)+m*U+dV/dy=0

i say (dU/dy)=W

so the precedent equations become

(dU/dy)=W

(dV/dy)=[(1-m)*(y/2)]*W-m*U

(dW/dy)=[V-(1-m)*(y/2)*U]*W+m*U^2-m

then on can write

(dU/dy)=f1(x[j],y[j])

(dV/dy)=f2(x[j],y[j])

(dW/dy)=f3(x[j],y[j])

right?

now with rk 4 i have:

U[j+1] = U[j] + (1/6)*(k1u + 2*k2u + 2*k3u + k4u);
where
k1u = h*f1(x[j],y[j]);
k2u = h*f1(x[j]+0.5*h,y[j] + 0.5*k1u);
k3u = h*f1(x[j] + 0.5*h,y[j] + 0.5*k2u);
k4u = h*f1(x[j] + h, y[j] + k3u);

V[j+1] = V[j] + (1/6)*(k1v + 2*k2v + 2*k3v + k4v);
where
k1v = h*f2(x[j],y[j]);
k2v = h*f2(x[j]+0.5*h,y[j] + 0.5*k1v);
k3v = h*f2(x[j] + 0.5*h,y[j] + 0.5*k2v);
k4v = h*f2(x[j] + h, y[j] + k3v);

W[j+1] = W[j] + (1/6)*(k1w + 2*k2w + 2*k3w + k4w);
where
k1w = h*f3(x[j],y[j]);
k2w = h*f3(x[j]+0.5*h,y[j] + 0.5*k1w);
k3w = h*f3(x[j] + 0.5*h,y[j] + 0.5*k2w);
k4w = h*f3(x[j] + h, y[j] + k3w);

so for my exercise i need to determine 12 value (k1u k2u k3u k4u k1v k2v k3v k4v k1w k2w k3w k4w) to determine U[j+1] V[j+1] and W[j+1]

right?

thanks

Do you know how to extend your system so that you use matrices instead of normal 1D derivatives?

no

Basically instead of having an f, you have a matrix A and a vector of f's (column vector) where you calculate A*v where v is the vector of f's.

The A matrix is the derivative matrix that contains all the derivative information for each f.