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