Clear Filters
Clear Filters

Function error in ODEs RK4 method, solving 6 unknowns

4 views (last 30 days)
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.

Accepted Answer

Alan Stevens
Alan Stevens on 8 Mar 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 Comments
VBBV
VBBV on 8 Mar 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)
Common
Common on 8 Mar 2024
Thanks, It works. Im a beginner, pls dont mind

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!