Function error in ODEs RK4 method, solving 6 unknowns
3 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Common
am 8 Mär. 2024
Kommentiert: Common
am 8 Mär. 2024
If I run this code, I get the following error.
Unrecognized function or variable 'func'.
Error in odearguments (line 92)
f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.
Error in ode45 (line 104)
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);
% Solves a system of IVP's using ode45
% Term Project (Dimensionless)
clear, clc
global CDin psi theta1 theta2 a b
fileID = fopen('Project_ode45.txt','w');
% Define Parameters for the problem
CDin = 4; %mol/L of Co2 inlet
psi = 0.1809; %psi=v/(V*k2)
theta1 = 0.0879*CDin^2; %(k1/k2)*Cdin^2 with
theta2 = (3.587/10000)*CDin; %k3/k2
a = 3; %N/C need to specify
b = 0.6; %H/C need to specify
y0 = [a;0;0;1;0;b]; %only co2 is 1, we can change a and b
tspan = [0 2];
[T,Y] = ode45(@func,tspan,y0);
dummy = [T,Y]; % Store output into a dummy matrix28
plot(dummy(:,1),dummy(:,2),dummy(:,1),dummy(:,3),dummy(:,1),dummy(:,4),dummy(:,1 ),dummy(:,5),dummy(:,1),dummy(:,6),dummy(:,1),dummy(:,7)) % Plot the solution
fprintf( ' Time y1 y2 y3 y4 y5 y6\n'); % Print output headings
fprintf(fileID,' Time y1 y2 y3 y4 y5 y6\r\n'); % Print output headings
fprintf( '%12.6f%12.6f%12.6f%12.6f%12.6f%12.6f%12.6f\n',dummy');
% Print to screen
fprintf(fileID,'%12.6f%12.6f%12.6f%12.6f%12.6f%12.6f%12.6f\r\n',dummy');
% Print to file
function f = func(t,y)
% function func to input the RHS of IVP's
global psi theta1 theta2 a b
f = zeros(length(y),1)+t*0; % Reset RHS functions
f(1) = psi*(a-y(1))-2*theta1*y(1)^2*y(4)+theta2*y(5)^2;
f(2) = theta2*y(5)^2-psi*y(2);
f(3) = theta1*y(1)^2*y(4)-y(3)-psi*y(3);
f(4) = psi*(1-y(4))-theta1*y(1)^2*y(4);
f(5) = y(3)-2*theta2*y(5)^2-psi*y(5);
f(6) = psi*(b-y(6))+y(3);
end
The following equations are need to be solved. y1 to y6 are the variables that are need to be find.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1638181/image.jpeg)
0 Kommentare
Akzeptierte Antwort
Alan Stevens
am 8 Mär. 2024
It works for me (I've commented out the fprint statements in the code below - but it works with them in!):
global CDin psi theta1 theta2 a b
% fileID = fopen('Project_ode45.txt','w');
% Define Parameters for the problem
CDin = 4; %mol/L of Co2 inlet
psi = 0.1809; %psi=v/(V*k2)
theta1 = 0.0879*CDin^2; %(k1/k2)*Cdin^2 with
theta2 = (3.587/10000)*CDin; %k3/k2
a = 3; %N/C need to specify
b = 0.6; %H/C need to specify
y0 = [a;0;0;1;0;b]; %only co2 is 1, we can change a and b
tspan = [0 2];
[T,Y] = ode45(@func,tspan,y0);
dummy = [T,Y]; % Store output into a dummy matrix28
plot(dummy(:,1),dummy(:,2),dummy(:,1),dummy(:,3),dummy(:,1),dummy(:,4),dummy(:,1 ),dummy(:,5),dummy(:,1),dummy(:,6),dummy(:,1),dummy(:,7)) % Plot the solution
% fprintf( ' Time y1 y2 y3 y4 y5 y6\n'); % Print output headings
% fprintf(fileID,' Time y1 y2 y3 y4 y5 y6\r\n'); % Print output headings
% fprintf( '%12.6f%12.6f%12.6f%12.6f%12.6f%12.6f%12.6f\n',dummy');
% % Print to screen
% fprintf(fileID,'%12.6f%12.6f%12.6f%12.6f%12.6f%12.6f%12.6f\r\n',dummy');
% Print to file
function f = func(t,y)
% function func to input the RHS of IVP's
global psi theta1 theta2 a b
f = zeros(length(y),1)+t*0; % Reset RHS functions
f(1) = psi*(a-y(1))-2*theta1*y(1)^2*y(4)+theta2*y(5)^2;
f(2) = theta2*y(5)^2-psi*y(2);
f(3) = theta1*y(1)^2*y(4)-y(3)-psi*y(3);
f(4) = psi*(1-y(4))-theta1*y(1)^2*y(4);
f(5) = y(3)-2*theta2*y(5)^2-psi*y(5);
f(6) = psi*(b-y(6))+y(3);
end
3 Kommentare
VBBV
am 8 Mär. 2024
Don't run the code in command window. Copy the code to a script file and save it using a valid filename, then run it from editor window using the Run button (green color)
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Ordinary Differential Equations finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!