why do I receive ode15s error?

Here is my code:
clc
clear
close all
%Data
k = 247; % W/mK
rho = 2710; % kg/m^3
Cp = 900; % J/kgK
h = 120; % W/m^2K
D = 0.005; % m
To = 298.15; % K
T_infinity = 298.15; % K
Tb = 423.15; % K
Nx = 100; % Step
xo = 0; % m
xL = 0.1; % m
% Step and Increment
x = linspace(xo, xL, Nx);
dx = x(2) - x(1); % delta-x
t = linspace(0,3600,100); % time till 3600 s (1 hour)
% Simplify Calculation
a = k/(rho*Cp);
b = 4*h./(D*rho*Cp);
c = 2*dx*h./k;
% Initial Condition
IC = To.*ones(Nx,1);
% Solver ODE15s
[t,T] = ode15s (@f,t,IC,[],xo,xL,D,Tb,T_infinity,To,Nx,a,b,c,dx,x);
% Recalculation
T(:,1) = Tb; %BC1
T(:,end) = (c.*T_infinity + 4*T(:,end-1)-T(:,end-2))./(3+c); %BC2
% Visualization
imagesc(x,t,T)
colormap jet
colorbar
grid on
title ('Temperature Profile')
xlabel ('length of the fin')
ylabel ('Time in seconds')
% Function
function dTdt = f(t,T,xo,xL,D,Tb,T_infinity,To,Nx,a,b,c,dx,x)
dTdt = zeros(Nx,1);
T(1) = Tb; %BC1
T(end)= (c.*T_infinity + 4*T(end-1) - T(end-2))./(3+c); %BC2
% Interior
for i = 2 : Nx-1
d2Tdx2(i) = (T(i+1)- 2*T(i) + T(i-1))./dx.^2;
dTdt = (a.*d2Tdx2(i) - b.*(T(i)-T_infinity));
end
end

 採用された回答

Torsten
Torsten 2022 年 10 月 25 日

1 投票

clc
clear
close all
%Data
k = 247; % W/mK
rho = 2710; % kg/m^3
Cp = 900; % J/kgK
h = 120; % W/m^2K
D = 0.005; % m
To = 298.15; % K
T_infinity = 298.15; % K
Tb = 423.15; % K
Nx = 100; % Step
xo = 0; % m
xL = 0.1; % m
% Step and Increment
x = linspace(xo, xL, Nx);
dx = x(2) - x(1); % delta-x
t = linspace(0,3600,100); % time till 3600 s (1 hour)
% Simplify Calculation
a = k/(rho*Cp);
b = 4*h./(D*rho*Cp);
c = 2*dx*h./k;
% Initial Condition
IC = To.*ones(Nx,1);
% Solver ODE15s
[t,T] = ode15s (@(t,T)f(t,T,xo,xL,D,Tb,T_infinity,To,Nx,a,b,c,dx,x),t,IC);
% Recalculation
T(:,1) = Tb; %BC1
T(:,end) = (c.*T_infinity + 4*T(:,end-1)-T(:,end-2))./(3+c); %BC2
% Visualization
imagesc(x,t,T)
colormap jet
colorbar
grid on
title ('Temperature Profile')
xlabel ('length of the fin')
ylabel ('Time in seconds')
% Function
function dTdt = f(t,T,xo,xL,D,Tb,T_infinity,To,Nx,a,b,c,dx,x)
dTdt = zeros(Nx,1);
T(1) = Tb; %BC1
T(end)= (c.*T_infinity + 4*T(end-1) - T(end-2))./(3+c); %BC2
% Interior
for i = 2 : Nx-1
d2Tdx2(i) = (T(i+1)- 2*T(i) + T(i-1))./dx.^2;
dTdt(i) = (a.*d2Tdx2(i) - b.*(T(i)-T_infinity));
end
end

1 件のコメント

RR
RR 2022 年 10 月 25 日
Thank you so much!

サインインしてコメントする。

その他の回答 (0 件)

カテゴリ

製品

リリース

R2022b

タグ

質問済み:

RR
2022 年 10 月 25 日

コメント済み:

RR
2022 年 10 月 25 日

Community Treasure Hunt

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

Start Hunting!

Translated by