I am trying to code a solution to blasius eq using Runge kutta 4, help please.

4 views (last 30 days)
clear all;
clc;
% 3 First order ODE´S from Blasius Eq
% dF/deta = G
% dG/deta = H
% dH/deta = -0.5*F*H
fF=@(eta,G) G;
fG=@(eta,H) H;
fH=@(eta,F,H) -0.5*F*H;
%initial conditions
F0 = 0;
G0 = 0;
H0 = 0; %Inital Guess for H0
% Step size and Eta max
h=0.0001;
eta=10;
N=ceil(eta/h);
%Update loop
for i=1:N
eta(i+1)=eta(i)+h;
% Runge-Kutta 4
k1F=fF(eta(i) ,F(i) ,G(i) ,H(i));
k1G=fG(eta(i) ,F(i) ,G(i) ,H(i));
k1H=fH(eta(i) ,F(i) ,G(i) ,H(i));
k2F=fF(eta(i)+h/2,F(i)+h/2*k1F,G(i)+h/2*k1G,H(i)+h/2*k1H);
k2G=fG(eta(i)+h/2,F(i)+h/2*k1F,G(i)+h/2*k1G,H(i)+h/2*k1H);
k2H=fH(eta(i)+h/2,F(i)+h/2*k1F,G(i)+h/2*k1G,H(i)+h/2*k1H);
k3F=fF(eta(i)+h/2,F(i)+h/2*k2F,G(i)+h/2*k2G,H(i)+h/2*k2H);
k3G=fG(eta(i)+h/2,F(i)+h/2*k2F,G(i)+h/2*k2G,H(i)+h/2*k2H);
k3H=fH(eta(i)+h/2,F(i)+h/2*k2F,G(i)+h/2*k2G,H(i)+h/2*k2H);
k4F=fF(eta(i)+h ,F(i)+h *k3F,G(i)+h *k3G,H(i)+h *k3H);
k4G=fG(eta(i)+h ,F(i)+h *k3F,G(i)+h *k3G,H(i)+h *k3H);
k4H=fH(eta(i)+h ,F(i)+h *k3F,G(i)+h *k3G,H(i)+h *k3H);
F(i+1)=F(i)+(h/6)*(k1F + 2*k2F + 2*k3F + k4F);
G(i+1)=G(i)+(h/6)*(k1G + 2*k2G + 2*k1G + k4G);
H(i+1)=H(i)+(h/6)*(k1G + 2*k2G + 2*k1G + k4G);
end
%Plot solution
figure(1); clf(1)
plot(eta,G)

Answers (1)

Torsten
Torsten on 3 Dec 2022
clear all;
clc;
% 3 First order ODE´S from Blasius Eq
% dF/deta = G
% dG/deta = H
% dH/deta = -0.5*F*H
fF=@(eta,F,G,H) G;
fG=@(eta,F,G,H) H;
fH=@(eta,F,G,H) -0.5*F*H;
%initial conditions
F0 = 0;
G0 = 0;
H0 = 0; %Inital Guess for H0
F(1) = F0;
G(1) = G0;
H(1) = H0;
% Step size and Eta max
h=0.0001;
eta=10;
N=ceil(eta/h);
%Update loop
for i=1:N
eta(i+1)=eta(i)+h;
% Runge-Kutta 4
k1F=fF(eta(i) ,F(i) ,G(i) ,H(i));
k1G=fG(eta(i) ,F(i) ,G(i) ,H(i));
k1H=fH(eta(i) ,F(i) ,G(i) ,H(i));
k2F=fF(eta(i)+h/2,F(i)+h/2*k1F,G(i)+h/2*k1G,H(i)+h/2*k1H);
k2G=fG(eta(i)+h/2,F(i)+h/2*k1F,G(i)+h/2*k1G,H(i)+h/2*k1H);
k2H=fH(eta(i)+h/2,F(i)+h/2*k1F,G(i)+h/2*k1G,H(i)+h/2*k1H);
k3F=fF(eta(i)+h/2,F(i)+h/2*k2F,G(i)+h/2*k2G,H(i)+h/2*k2H);
k3G=fG(eta(i)+h/2,F(i)+h/2*k2F,G(i)+h/2*k2G,H(i)+h/2*k2H);
k3H=fH(eta(i)+h/2,F(i)+h/2*k2F,G(i)+h/2*k2G,H(i)+h/2*k2H);
k4F=fF(eta(i)+h ,F(i)+h *k3F,G(i)+h *k3G,H(i)+h *k3H);
k4G=fG(eta(i)+h ,F(i)+h *k3F,G(i)+h *k3G,H(i)+h *k3H);
k4H=fH(eta(i)+h ,F(i)+h *k3F,G(i)+h *k3G,H(i)+h *k3H);
F(i+1)=F(i)+(h/6)*(k1F + 2*k2F + 2*k3F + k4F);
G(i+1)=G(i)+(h/6)*(k1G + 2*k2G + 2*k3G + k4G);
H(i+1)=H(i)+(h/6)*(k1H + 2*k2H + 2*k3H + k4H);
end
%Plot solution
figure(1); clf(1)
plot(eta,G)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!