问题描述：

反应扩散方程Matlab编程

>> function fd1d_predator_prey ( )

% FD1D_PREDATOR_PREY.m one-dimensional finite-difference method for Scheme 2

% applied to the predator-prey system with Kinetics 1.

%

% Author:

%

% Marcus Garvie

%

%

% User inputs of parameters

%

alpha = input('Enter parameter alpha ');

beta = input('Enter parameter beta ');

gamma = input('Enter parameter gamma ');

delta = input('Enter parameter delta ');

a = input('Enter a in [a,b] ');

b = input('Enter b in [a,b] ');

h = input('Enter space-step h ');

T = input('Enter maximum time T ');

delt = input('Enter time-step Delta t ');

%

% User inputs of initial data

%

u0 = input('Enter initial data function u0(x) ','s'); % see notes

v0 = input('Enter initial data function v0(x) ','s'); % in text

%

% Calculate some constants

%

mu=delt/(h^2);

J=round((b-a)/h);

n = J+1; % no.of nodes (d.f.) for each dependent variable

N=round(T/delt);

%

% Initialization

%

u=zeros(n,1); v=zeros(n,1); F=zeros(n,1); G=zeros(n,1);

y1=zeros(n,1); y2=zeros(n,1); z1=zeros(n,1); z2=zeros(n,1);

B1=sparse(n,n); B2=sparse(n,n); L=sparse(n,n); Lower1=sparse(n,n);

Lower2=sparse(n,n); U1=sparse(n,n); U2=sparse(n,n);

%

% Assign initial data

%

indexI=[1:n]';

x=(indexI-1)*h+a; % vector of x values on grid

u = eval(u0).*ones(n,1); v = eval(v0).*ones(n,1);

%

% Construct matrix L (without 1/h^2 factor)

%

L=sparse(1,2,-2,n,n);

L=L+sparse(n,n-1,-2,n,n);

L=L+sparse(2:n-1,3:n,-1,n,n);

L=L+sparse(2:n-1,1:n-2,-1,n,n);

L=L+sparse(1:n,1:n,2,n,n);

%

% Construct matrices B1 & B2

%

B1=sparse(1:n,1:n,1,n,n) + mu*L;

B2=sparse(1:n,1:n,1,n,n) + delta*mu*L;

%

% Perform the LU factorisation of B1 and B2

%

[Lower1,Upper1]=lu(B1);

[Lower2,Upper2]=lu(B2);

%

% Time-stepping procedure

%

for nt=1:N

% Evaluate modified functional response

hhat = u./(alpha + abs(u));

% Update right-hand-side of linear system

F = u - u.*abs(u) - v.*hhat;

G = beta*v.*hhat - gamma*v;

y1 = u + delt*F;

y2 = v + delt*G;

% Forward substitution to solve Lower1*z1=y1 for z1

z1 = Lower1\y1;

% Back-substitution to solve Upper1*u=z1 for u

u = Upper1\z1;

% Forward substitution to solve Lower2*z2=y2 for z2

z2 = Lower2\y2;

% Back-substitution to solve Upper2*v=z2 for v

v = Upper2\z2;

end

%

% Plot solution at time level T=N*delt

%

plot(x,u,'k'); hold on; plot(x,v,'k-.')

return

end

function fd1d_predator_prey ( )

|

Error:Function definitions are not permitted at the prompt or in scripts.

>>

1个回答

分类：

综合

2014-10-06

问题解答：

我来补答

1

将代码保存为M文件

2

将下面的代码复制到命令窗口,按回车键,运行即可

fd1d_predator_prey

展开全文阅读