Error from ode15s
2 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Received this error: Error in ode15s (line 153)
odearguments(odeIsFuncHandle, odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);
Below is my code:
clc
psi = 0.63;
psip = 0.27;
L = 150e-6;
Ts = 60;
k = 0.00158;
q0 = 4.62e3;
C0 =1.03;
gama = 5021.8;
Jw = 8.5e-6;
a = psi*L/Jw/Ts;
b = (1-psi-psip)*k*L/Jw;
e = q0./C0;
d = k*Ts*C0./q0;
%TIME SPAN
tt = 0:100:2400;
t = tt/Ts;
%STEP LENGTH
Nz = 51;
zz= linspace(0, L, Nz);
z = zz./L;
dz = z(2) -z(1);
%IC
IC_A = zeros(1, Nz);
IC_B = zeros(1, Nz);
IC = [IC_A, IC_B];
%Solve the pde using ode15s
[t, c] = ode15s(@f, t, IC, [], Nz, C0, a, b, e, d, dz, gama);
% Extract solutions
C = c(:, 1:Nz);
g = c(:, Nz+1 : 2*Nz);
% Reinput BC
C(:, 1) = 1;
C(:, end) = (4*C(:, Nz-1)-C(:, Nz-2))./3;
%plot
figure(1)
plot(t, C(:, end), '--')
ylim([0 1.2])
function dcqdt = f(t, c, Nz, C0, a, b, e, d, dz, gama)
dcqdt = zeros(length(c), 1);
dcdt = zeros(Nz, 1);
dqdt = zeros(Nz, 1);
%define value
C = c(1:Nz);
g = c(Nz+1 : 2*Nz);
%BC
C(1) =1;
C(end) = (4*C(Nz-1)-C(Nz-2))./3;
%interior for-loop
for i=2:Nz-1
dcdz(i) = C(i+1) - C(i-1)./2./dz; % Using FDA
dcdt(i) = (-1./a).*(dcdz(i) + b.*(gama.*C(i)-q(i).*e));
dqdt(i) = d.*(gama.*C(i)-q(i).*e);
end
dcqdt = [dcdt; dqdt];
end
0 commentaires
Réponses (1)
Torsten
le 27 Juil 2022
Modifié(e) : Torsten
le 27 Juil 2022
psi = 0.63;
psip = 0.27;
L = 150e-6;
Ts = 60;
k = 0.00158;
q0 = 4.62e3;
C0 =1.03;
gama = 5021.8;
Jw = 8.5e-6;
a = psi*L/Jw/Ts;
b = (1-psi-psip)*k*L/Jw;
e = q0./C0;
d = k*Ts*C0./q0;
%TIME SPAN
tt = 0:100:24000;
t = tt/Ts;
%STEP LENGTH
Nz = 51;
zz= linspace(0, L, Nz);
z = zz./L;
dz = z(2) -z(1);
%IC
IC_A = zeros(1, Nz);
IC_B = zeros(1, Nz);
IC = [IC_A, IC_B];
%Solve the pde using ode15s
[t, c] = ode15s(@(t,y)f(t, y, Nz, C0, a, b, e, d, dz, gama), t, IC);
% Extract solutions
C = c(:, 1:Nz);
g = c(:, Nz+1 : 2*Nz);
% Reinput BC
C(:, 1) = 1;
C(:, end) = (4*C(:, Nz-1)-C(:, Nz-2))./3;
%plot
figure(1)
plot(t, C(:, end), '--')
%ylim([0 1.2])
function dcqdt = f(t, c, Nz, C0, a, b, e, d, dz, gama)
dcqdt = zeros(length(c), 1);
dcdt = zeros(Nz, 1);
dqdt = zeros(Nz, 1);
%define value
C = c(1:Nz);
q = c(Nz+1 : 2*Nz);
%BC
C(1) =1;
C(end) = (4*C(Nz-1)-C(Nz-2))./3;
%interior for-loop
for i=2:Nz-1
dcdz(i) = (C(i+1) - C(i-1))./2./dz; % Using FDA
dcdt(i) = (-1./a).*(dcdz(i) + b.*(gama.*C(i)-q(i).*e));
dqdt(i) = d.*(gama.*C(i)-q(i).*e);
end
dcqdt = [dcdt; dqdt];
end
Voir également
Catégories
En savoir plus sur Geometry and Mesh dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!