Solving system of ODEs with ode45
1 view (last 30 days)
Show older comments
i have the following equetions system
and and the maching valuse
i tried to solve the system with the following code but i dont get the right resultes which i know:
ephsilon values [1.94 2.65 3.42] for the x1/h vec [7.5 9.5 11]
k values [0.2680 0.3585 0.444] for the same x1/h vec.
the code i wrote
Y=zeros(1,2);
Uc=12.4; %m/s
S=46.8;% 1/s
k_75=0.268 ; %m2/s2;
e_75=1.94; %m2/s3
Cu = 0.09;
C1= 1.44;
C2=1.92;
x1_h=7.5;
h = 0.3048;%m
xspan=[7.5 9.5 11];
Y0=[k_75 e_75 ];
f = @(t,Y) [(Cu*((Y(1))^2)/(Y(2))*(S^2)-Y(2))/(Uc/h); (Cu*C1*(Y(1))*(S^2)-C2*((Y(2))^2)/(Y(1)))/(Uc/h)];
[x0,val] = ode45(f , xspan, Y0);
figure (1)
plot(x0, val(:,1))
figure (2)
plot(x0, val(:,2))
what wrong with it ? thanks for the help
1 Comment
Answers (2)
Sam Chak
on 24 Dec 2022
Edited: Sam Chak
on 24 Dec 2022
Hi @shir levi
Edit: Previously I've got the same results like you did. However, after thinking about it, those 2 sets of 3 values for k and ϵ are probably initial conditions. So, I modified the code to include a for loop.
Also, in the 3rd figure, the system is shown to be inherently unstable, and there is no equilibrium point.
% parameters
Uc = 12.4; % m/s
S = 46.8; % 1/s
k = [0.2680 0.3585 0.444]; % m2/s2 ic for k
epsil = [1.9400 2.6500 3.420]; % m2/s3 ic for epsil
Cu = 0.09;
C1 = 1.44;
C2 = 1.92;
x1_h = [7.5 9.5 11];
h = 0.3048; % m
tstart = x1_h*h/Uc;
tfinal = 2*tstart;
for j = 1:length(x1_h)
% setting conditions
f = @(t, Y) [Cu*(Y(1)^2)/Y(2)*S^2 - Y(2); Cu*C1*Y(1)*S^2 - C2*(Y(2)^2)/Y(1)];
tspan = [tstart(j) tfinal(j)];
Y0 = [k(j) epsil(j)];
[t, Y] = ode45(f, tspan, Y0);
% plotting results
figure(1)
plot(t*Uc/h, Y(:,1)), hold on, grid on,
xlabel('x_{1}/h'), ylabel('k'), title('k responses')
figure(2)
plot(t*Uc/h, Y(:,2)), hold on, grid on,
xlabel('x_{1}/h'), ylabel('\epsilon'), title('\epsilon responses')
end
hold off
[x, y] = meshgrid(-11:2:11);
u = Cu*(x.^2)./y*S^2 - y;
v = Cu*C1*x*S^2 - C2*(y.^2)./x;
figure(3)
l = streamslice(x, y, u, v);
axis tight
xlabel('k')
ylabel('\epsilon')
0 Comments
shir levi
on 24 Dec 2022
1 Comment
Sam Chak
on 25 Dec 2022
Hi @shir levi
But you used these values and in the initial condition.
So, I'm confused.
Perhaps, if you tell us more info about the 2 sets of 3 values, it will help clarifying the matter.
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!