Non-uniform grid

3 次查看（过去 30 天）

颯太小濱 2022-1-18

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1631065-non-uniform-grid

评论： John D'Errico 2022-3-9

I have a question about unequal spacing.

In the code below, I would like to either increase the number of N to avoid running out of memory, or increase the number of grids only around i=2N to prevent the solution from diverging.

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(321221,1,N)];
options = odeset('Vectorized','on');
[t,y] = ode15s(@f,tspan,y0,options);
u= y(:,1:2:end);
x = (1:N)/(N+1);
p=y(:,2:2:end);
assignin('base','s',u(end,:))
assignin('base','p',p(end,:))
    function dydt = f(~,y)
        a=0.0002/N;
        b=3.15*10^-4;
        c=2.4*10^-5;
        d=300000;
        e=1;
        f=174284;
        k=6.4*10^-14; 
        l=5;
       
        dydt = zeros(2*N,size(y,2));      % preallocate dy/dt
          
        i = 1;
        dydt(i,:) = 1/2/a^2/b*((k*(y(i,:).^l+y(i+2,:).^l)).*((y(i+3,:)-y(i+1,:))+f*(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*(2* y(i,:).^l)).*((y(i+1,:)-d)+f*(1.417*(y(i,:)-e)-2.12*((1-e).^2-(1-y(i,:)).^2)+1.263*((1-e).^3-(1-y(i,:)).^3))).*y(i,:));
        dydt(i+1,:) =1/2/a^2/c*((k*((1-y(i,:)).^l+(1-y(i+2,:)).^l).*(y(i+3,:)-y(i+1,:))).*y(i+3,:)-(2*k*((1-y(i,:)).^l).*(y(i+1,:)-d)).*y(i+1,:)); 
        % 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,:).^l+y(i+2,:).^l)).*((y(i+3,:)-y(i+1,:))+f*(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,:).^l+y(i,:).^l)).*((y(i+1,:)-y(i-1,:))+f*(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,:)).^l+(1-y(i+2,:)).^l).*(y(i+3,:)-y(i+1,:))).*y(i+3,:)-(k*((1-y(i-2,:)).^l+(1-y(i,:)).^l).*(y(i+1,:)-y(i-1,:))).*y(i+1,:));
        
        i = 2*N-1;
        dydt(i,:) = 1/2/a^2/b*(-9.93*10^-6*a*b*2-(k*(y(i-2,:).^l+y(i,:).^l)).*((y(i+1,:)-y(i-1,:))+f*(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*(804./y(i,:) *a*c*2-(k*((1-y(i-2,:)).^l+(1-y(i,:)).^l)).*(y(i+1,:)-y(i-1,:)).*y(i+1,:));
        
    end 
end