Can someone help me to correct the code ​​for this problem using ode45 solver?

1 Ansicht (letzte 30 Tage)
VAIDEHI
VAIDEHI am 9 Mär. 2023
Kommentiert: Alex Sha am 12 Mär. 2023
Question is:
Governing equations:
Boundary conditions:
f'(10)= 0, (assume )
As tried here, I want to plot graphs for velocity profile using parameters like Pr, beta, and N.
I tried giving 'for' loop, but it failed.
How to enter the code command to plot a graph for velocity profile with different beta parameter?
Trial code
function dy=vj_Casson()
clc; clear all;
betava=[0.5 1 2 5];
for i=1:4
beta=betava(i);
tspan = [0 10];
S = 0.5;
y0 = [S 0 0 1 1];
% y0 = [S 1 0 1 0];
[eta,y] = ode45(@fun,tspan,y0);
plot(eta, y(:,1));
xlabel('\bf\eta','FontSize',20,'FontWeight','bold');
ylabel('f(\eta)','FontSize',10,'FontWeight','bold');
legend('\beta=0.5','\beta=1','\beta=2','\beta=5')
hold on
end
end
function dy = fun(eta,y)
dy = zeros(5,1);
Pr=0.3; N=0.1; beta=0.5;
dy(1) = y(2);
dy(2)=y(3);
dy(3)=(((2*y(2))^(2))-y(1)*y(3))./((1+(1./beta)));
dy(4) =y(5);
dy(5)=-(((3*Pr)*(y(1)*y(5)-y(2)*y(4)))/(4*N+3));
end
  1 Kommentar
Alex Sha
Alex Sha am 12 Mär. 2023
The results below are what you want?
1: f' = df/dt = f'
2: f'' = df'/dt = f''
3: theta' = dtheta/dt = theta'
4: f''' = df''/dt = (2*(f')^2-f*f'')/(1+1/0.5)
5: theta'' = dtheta'/dt = 0.3*(f'*theta-f*theta')/(1+4/3*0.1)
Objective Function: 1.41854343905794E-26
Boundary Values Estimated:
f''(t=0): -0.887446621710334
theta'(t=0): -0.36495059823582

Melden Sie sich an, um zu kommentieren.

Antworten (1)

Torsten
Torsten am 9 Mär. 2023
This is a boundary value problem, not an initial value problem since conditions on f and theta are given on both ends of the integration interval (eta = 0 and eta = 10). You have to use bvp4c or bvp5c instead of ode45 to solve.
  3 Kommentare
Torsten
Torsten am 10 Mär. 2023
Bearbeitet: Torsten am 10 Mär. 2023
You cannot use conditions at eta = 10 if you use ode45. You must assume two further conditions at eta = 0 to make ode45 work and try to adjust these conditions in several runs such that you arrive at your two conditions at eta = 10. Look up "shooting method" for more details.
I can assure you: the simpler way for you to go is to use bvp4c or bvp5c.

Melden Sie sich an, um zu kommentieren.

Community Treasure Hunt

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

Start Hunting!

Translated by