Setting initial condition as an inverse function in ODE15s wont plot. Am I missing something?

Messing around with the diffusion equation discretized to du(n)/dt=K(u(n+1)-2u(n)+u(n-1) and changing the value of K but when I try to set K= 1/u(n) a long error message occurs with: Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN. Warning: Failure at t=0.000000e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (7.905050e-323) at time t. Code Provided:
function [t,u]=discretized
tspan = [0 100];
n=100;
y0=zeros(n,1);
for i=1:n
y0(n/2)=1;
end
[t,u]=ode15s(@f,tspan,y0);
for i=1:length(t)
plot(u(i,:))
title('Discretization of Diffusion Equation');
xlabel('Cell');
ylabel('Concentration');
pause(1);
end
end
function dudt=f(t,y)
n=100;
dudt=zeros(n,1);
K=1./(y(1));
dudt(1)=K*(y(2)-y(1));
for i=2:n-1
K=1./(y(i));
dudt(i)=K*(y(i+1)-2*(y(i))+y(i-1));
end
dudt(n)=K*(y(n-1)-y(n));
end
Any ideas on what the problem is?

6 件のコメント

Torsten
Torsten 2018 年 2 月 26 日
Since y0=0 almost everywhere, you divide by 0 almost everywhere.
Best wishes
Torsten.
Philip Mitchell
Philip Mitchell 2018 年 2 月 26 日
Thats the initial condition for only the middle cell to have any concentration. I set the middle cell to 1 and it should spread out. Thanks
Torsten
Torsten 2018 年 2 月 26 日
Yes, but in the calculation of K, you divide by the concentrations of the cells, and these are 0.
Philip Mitchell
Philip Mitchell 2018 年 2 月 26 日
Ah nice spot! So I would remove the line dudt=zeros(n,1); or alter it to something else? Thanks again!
Torsten
Torsten 2018 年 2 月 26 日
No. The initial concentration y0 can't be 0 if you set K=1/y.
Philip Mitchell
Philip Mitchell 2018 年 2 月 26 日
Set the matrix to a n*1 matrix of values close to zero and this now works I think thanks!

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

回答 (0 件)

カテゴリ

ヘルプ センター および File ExchangeMathematics についてさらに検索

製品

質問済み:

2018 年 2 月 26 日

コメント済み:

2018 年 2 月 26 日

Community Treasure Hunt

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

Start Hunting!

Translated by