memory shortage: Solve Stiff ODEs
2 views (last 30 days)
Show older comments
Dear matlab users.
I have memory shortage problem .
I use Solve Stiff ODEs https://jp.mathworks.com/help/matlab/math/solve-stiff-odes.html?lang=en
function brussode(N)
if nargin<1
N = 100;
end
tspan = [0; 10];
y0 = [repmat(0.9,1,N); repmat(100200,1,N)];
options = odeset('Vectorized','on','JPattern',jpattern(N));
[t,y] = ode15s(@f,tspan,y0,options);
u = y(:,1:2:end);
x = (1:N)/(N+1);
p=y(:,2:2:end);
figure;
plot(x,u(end,:))
figure;
plot(x,p(end,:))
function dydt = f(~,y)
a=0.0002/N;
b=3.55*10^-4;
c=2.39*10^-5;
k=3*10^-12;
dydt = zeros(2*N,size(y,2)); % preallocate dy/dt
i = 1;
dydt(i,:) = 1/2/a^2/b*((k*(y(i,:).^3+y(i+2,:).^3)).*((y(i+3,:)-y(i+1,:))+ +19206*(1.417*(y(i+2,:)-y(i,:))-2.12*((1-y(i,:)).^2-(1-y(i+2,:)).^2)+1.263*((1-y(i,:)).^3-(1-y(i+2,:)).^3))).*y(i,:)-1.29*10^-9*a*b*2);
dydt(i+1,:) =1/2/a^2/c*((k*((1-y(i,:)).^3+(1-y(i+2,:)).^3).*(y(i+3,:)-y(i+1,:))).*y(i+3,:)+0.1./y(i+1,:) *a*c*2);
% Evaluate the 2 components of the function at all interior grid points.
i = 3:2:2*N-3;
dydt(i,:) = 1/2/a^2/b*((k*(y(i,:).^3+y(i+2,:).^3)).*((y(i+3,:)-y(i+1,:))+19206*(1.417*(y(i+2,:)-y(i,:))-2.12*((1-y(i,:)).^2-(1-y(i+2,:)).^2)+1.263*((1-y(i,:)).^3-(1-y(i+2,:)).^3))).*y(i,:)-(k*(y(i-2,:).^3+y(i,:).^3)).*((y(i+1,:)-y(i-1,:))+19206*(1.417*(y(i,:)-y(i-2,:))-2.12*((1-y(i-2,:)).^2-(1-y(i,:)).^2)+1.263*((1-y(i-2,:)).^3-(1-y(i,:)).^3))).*y(i-2,:));
dydt(i+1,:) =1/2/a^2/c*((k*((1-y(i,:)).^3+(1-y(i+2,:)).^3).*(y(i+3,:)-y(i+1,:))).*y(i+3,:)-(k*((1-y(i-2,:)).^3+(1-y(i,:)).^3)).*(y(i+1,:)-y(i-1,:)).*y(i+1,:));
i = 2*N-1;
dydt(i,:) = 1/2/a^2/b*(1.29*10^-9*a*b*2-(k*(y(i-2,:).^3+y(i,:).^3)).*((y(i+1,:)-y(i-1,:))+19206*(1.417*(y(i,:)-y(i-2,:))-2.12*((1-y(i-2,:)).^2-(1-y(i,:)).^2)+1.263*((1-y(i-2,:)).^3-(1-y(i,:)).^3))).*y(i-2,:));
dydt(i+1,:) =1/2/a^2/c*(0.1./y(i+1,:) *a*c*2-(k*((1-y(i-2,:)).^3+(1-y(i,:)).^3)).*(y(i+1,:)-y(i-1,:)).*y(i+1,:));
;
end
end
function S = jpattern(N)
% Jacobian sparsity pattern
B = ones(2*N,5);
B(2:2:2*N,2) = zeros(N,1);
B(1:2:2*N-1,4) = zeros(N,1);
S = spdiags(B,-2:2,2*N,2*N);
end
Please let me know what I need to change to calculate under 300GB.
Thank you for your help.
0 Comments
Accepted Answer
Wan Ji
on 18 Aug 2021
你好!这个直接把
options = odeset('Vectorized','on','JPattern',jpattern(N));
改写成
options = odeset('Vectorized','on');
就可以了,因为你给出的ode方程的jpattern函数不能照搬别人的,所以干脆不用,这样计算也能快速得到想要的结果。
More Answers (0)
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!