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

6 次查看(过去 30 天)
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)

回答(2 个)

Torsten
Torsten 2022-12-3
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)

VBBV
VBBV 2024-9-9
@Guillermo, The anonymous functions, F ,G, H defined for the blasius flow need to applied in the same manner when RK4 method is implemented
clear all;
clc;
% 3 First order ODE´S from Blasius Eq
% dF/deta = G
% dG/deta = H
% dH/deta = -0.5*F*H
%initial conditions
F(1) = 0.01;
G(1) = 0.01;
H(1) = 0.1; %Inital Guess for H0
fF=@(eta,G) G;
fG=@(eta,H) H;
fH=@(eta,F,H) -0.5*F*H;
% 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),G(i));
k1G=fG(eta(i),H(i));
k1H=fH(eta(i),F(i),H(i));
k2F=fF(eta(i)+h/2,G(i)+h/2*k1G);
k2G=fG(eta(i)+h/2,H(i)+h/2*k1H);
k2H=fH(eta(i)+h/2,F(i)+h/2*k1F,H(i)+h/2*k1H);
k3F=fF(eta(i)+h/2,G(i)+h/2*k2G);
k3G=fG(eta(i)+h/2,H(i)+h/2*k2H);
k3H=fH(eta(i)+h/2,F(i)+h/2*k2F,H(i)+h/2*k2H);
k4F=fF(eta(i)+h,G(i)+h*k3G);
k4G=fG(eta(i)+h,H(i)+h*k3H);
k4H=fH(eta(i)+h,F(i)+h*k3F,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)*(k1H + 2*k2H + 2*k1H + k4H);
end
%Plot solution
hold on
subplot(311);plot(eta,F); subplot(312); plot(eta,G); subplot(313);plot(eta,H);

类别

Help CenterFile Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by