Piotr Krzyżanowski

Poniżej przestawiamy wybrane kody źródłowe programów omawianych w pierwszej części podręcznika Obliczenia inżynierskie i naukowe (PWN, 2011). Większość z nich jest dostępna na licencji Creative Commons BY-SA.

Listing 5.1 - invlaplace.m. [Pobierz]

%%
%% Listing 5.1 w książce "Obliczenia inżynierskie i naukowe", PWN 2011
%% (C) Piotr Krzyżanowski, 2011
%%
%% Możesz używać tego kodu zgodnie z licencją Creative Commons BY-SA, 
%% zob. http://creativecommons.org/licenses/by-sa/3.0/pl/
%%

% skrypt invlaplace.m
for i=1:N
	for j=1:i
		L(i,j) = j;
	end
	for j=i:N
		L(i,j) = i;
	end
end

Listing 5.2 - invlaplace2.m. [Pobierz]

%%
%% Listing 5.2 w książce "Obliczenia inżynierskie i naukowe", PWN 2011
%% (C) Piotr Krzyżanowski, 2011
%%
%% Możesz używać tego kodu zgodnie z licencją Creative Commons BY-SA, 
%% zob. http://creativecommons.org/licenses/by-sa/3.0/pl/
%%

% skrypt invlaplace.m - wersja zmodyfikowana
L = zeros(N,N);
for i=1:N
	for j=1:i
		L(i,j) = j;
	end
	for j=i:N
		L(i,j) = i;
	end
end

Listing 5.3 - invlaplace3.m. [Pobierz]

%%
%% Listing 5.2 w książce "Obliczenia inżynierskie i naukowe", PWN 2011
%% (C) Piotr Krzyżanowski, 2011
%%
%% Możesz używać tego kodu zgodnie z licencją Creative Commons BY-SA, 
%% zob. http://creativecommons.org/licenses/by-sa/3.0/pl/
%%

% skrypt invlaplace.m - wersja działająca na blokach macierzy
L = zeros(N,N);
for i=1:N
	L(i,i:N) = L(i:N,i) = i;
end

Listing 5.4 - invlap1d.m. [Pobierz]

%%
%% Listing 5.4 w książce "Obliczenia inżynierskie i naukowe", PWN 2011
%% (C) Piotr Krzyżanowski, 2011
%%
%% Możesz używać tego kodu zgodnie z licencją Creative Commons BY-SA, 
%% zob. http://creativecommons.org/licenses/by-sa/3.0/pl/
%%

function L = invlap1d(N)
% funkcja invlap1d - generuje macierz odwrotną do 1-wym laplasjanu 
% z warunkami Dirichleta z lewego i Neumanna z prawego krańca odcinka
%
% wywołanie: L = invlap1d(N)
%
% N - wymiar macierzy
% L - macierz odwrotna laplasjanu

L = zeros(N,N);
for i=1:N
	L(i,i:N) = L(i:N,i) = i;
end
end

Listing 5.5 - invlap1dsuper.m. [Pobierz]

%%
%% Listing 5.5 w książce "Obliczenia inżynierskie i naukowe", PWN 2011
%% (C) Piotr Krzyżanowski, 2011
%%
%% Możesz używać tego kodu zgodnie z licencją Creative Commons BY-SA, 
%% zob. http://creativecommons.org/licenses/by-sa/3.0/pl/
%%

function [L, norma] = invlap1d(N,p)
% Funkcja invlap1d - generuje macierz odwrotną do 1-wym laplasjanu 
% z warunkami Dirichleta z lewego i Neumanna z prawego krańca odcinka
%
% Wywołanie: 
%		[L, norma] = invlap1d(N, p)
% Parametry:
% 		N - rozmiar macierzy 
% 		p - jaka norma macierzy (opcjonalnie, domyślnie: 2)
% Wyniki:
% 		L - macierz odwrotna laplasjanu
% 		norma - (p)-ta norma macierzy L (opcjonalnie)

% obsługa parametrów wejściowych - wartości domyślne
if (nargin < 2)	
	p = 2;
	if (nargin < 1)
		error('Za mała liczba argumentów');
	end
end

% właściwa treść funkcji
L = zeros(N,N);
for i=1:N
	L(i,i:N) = L(i:N,i) = i;
end

% normę liczymy tylko wtedy, gdy użytkownik jej zażąda
if (nargout > 1)
	norma = norm(L,p);
end
end

Listing 5.7 - divdif.m. [Pobierz]

%%
%% Listing 5.7 w książce "Obliczenia inżynierskie i naukowe", PWN 2011
%% (C) Piotr Krzyżanowski, 2011
%%
%% Możesz używać tego kodu zgodnie z licencją Creative Commons BY-SA, 
%% zob. http://creativecommons.org/licenses/by-sa/3.0/pl/
%%

function df = divdif(f, x, h)
% 
% DF = divdif(F, X, H)
%
% Wyznacza przybliżenie F'(X) ze wzoru (F(X+H)-F(X-H)) / (2H)
% Jeśli nie podano H, przyjmuje się H = max(|X|,1)*sqrt(eps).

if nargin < 3
	h = max(abs(x),1)*sqrt(eps);
end
f1 = feval(f,x+h);
f2 = feval(f,x-h);
df = (f1-f2)./(2*h);
end

Listing 6.1 - funkcja2d.m. [Pobierz]

%%
%% Listing 6.1 w książce "Obliczenia inżynierskie i naukowe", PWN 2011
%% (C) Piotr Krzyżanowski, 2011
%%
%% Możesz używać tego kodu zgodnie z licencją Creative Commons BY-SA, 
%% zob. http://creativecommons.org/licenses/by-sa/3.0/pl/
%%

disp('Wizualizacja funkcji 2D'); 

% definicja funkcji do zwizualizowania: skorzystamy z funkcji anonimowej;
% dzięki użyciu operatorów tablicowych, argumenty x,y mogą być macierzami
ff = @(x,y) abs(sin(x-1)).^((y+1).^2).*abs(cos((y+1))).^((x-1).^2); 

% utworzenie siatki na dziedzinie funkcji
x = linspace(-3,3,30);  % węzły na osi OX
y = linspace(-4,4,40);  % węzły na osi OY
[X,Y] = meshgrid(x,y);  % siatka na płaszczyźnie

% wyznaczenie wartości ff na siatce
Z = ff(X,Y);

% wizualizacja:

contour(X,Y,Z);		% poziomice
pause(5);

mesh(X,Y,Z);		% powierzchnia
pause(5);

surf(X,Y,Z);		% powierzchnia inaczej
pause(5);

imagesc(X,Y,Z);		% mapa
axis('xy','square');
pause(5);

Listing 6.2 - stationary.m. [Pobierz]

%%
%% Listing 6.2 w książce "Obliczenia inżynierskie i naukowe", PWN 2011
%% (C) Piotr Krzyżanowski, 2011
%%
%% Możesz używać tego kodu zgodnie z licencją Creative Commons BY-SA, 
%% zob. http://creativecommons.org/licenses/by-sa/3.0/pl/
%%

clear;

delta = 9/7;
TOL = 1e-3; % jaką wartość $||F||^2$ uznajemy za "małą"?

% Definicja funkcji prawej strony przydatna do zabawy z polem wprost w skrypcie
% jest dopuszczalna tylko w Octave: MATLAB wymaga zdefiniowania tej funkcji 
% w osobnym pliku
function [f1, f2, f3] = pole(u1, u2, u3, delta)
f1 = 1 - u1 - delta*(u2.^2+u3.^2);
f2 = -u2 + delta*u1.*u2;
f3 = -u3 + delta*u1.*u3;
end

% szukamy punktów stacjonarnych w prostokącie $(u_2,u_3)\in [-A,A]\times [-A,A]$...
N = 30; A = 0.6;
u2 = linspace(-A,A,N);
u3 = linspace(-A,A,N);
[U2, U3] = meshgrid(u2,u3);

% ...i dla u1 = [0.7..0.8]
for u1 = linspace(0.7,1.2,50)
	[F1 F2 F3] = pole(u1, U2, U3, delta);
	G = double(F1.^2 + F2.^2 + F3.^2 < TOL);
	% czy znaleźliśmy choć jeden punkt, w którym F jest małe?
	if max(G(:)) > 0
		fprintf(stderr,'Mala norma dla u1 = %g',u1); % MATLAB nie zna stderr...
		[i, j] = find(G);
		fprintf(stderr,' i np. (u2,u3) = (%5.3g,%5.3g)\n', ...
					U2(i(1),j(1)), U3(i(1),j(1)));
		mesh(U2,U3,G); 
		title(['u1=', num2str(u1)]);
		pause(1);
	end
end

Listing 6.3 - spyc.m. [Pobierz]

%%
%% Listing 6.3 w książce "Obliczenia inżynierskie i naukowe", PWN 2011
%% (C) Piotr Krzyżanowski, 2011
%%
%% Możesz używać tego kodu zgodnie z licencją Creative Commons BY-SA, 
%% zob. http://creativecommons.org/licenses/by-sa/3.0/pl/
%%

function spyc(D)
if (nnz(D) < 10000) % wybór markera
	m ='*'; 
else
	m ='.';
end 
if nnz(D>0)
	spy(D>0, ['r',m]); % wizualizacja dodatnich elementów macierzy D 
	hold on;	
end
if nnz(D<0)
	spy(D<0, ['b',m]);	% wizualizacja ujemnych elementów macierzy D
end
hold off;
end

Listing 7.2 - solvegr.m. [Pobierz]

%%
%% Listing 7.2 w książce "Obliczenia inżynierskie i naukowe", PWN 2011
%% (C) Piotr Krzyżanowski, 2011
%%
%% Możesz używać tego kodu zgodnie z licencją Creative Commons BY-SA, 
%% zob. http://creativecommons.org/licenses/by-sa/3.0/pl/
%%

function x = solvegr(e, W, D, b, tol, maxit)
%
% Rozwiazuje (I - e*W*D^{-1})x = (1-e)b
% 	0 <= e <= 1, 
% 	W - symetryczna
% 	D - diagonalna, dod. okresl.

if nargin < 5
	tol = 1e-8;
end
if nargin < 6
	maxit = 100;
end

switch e
case 0
	x = b;
case 1
	opts.maxit = maxit; opts.tol = tol; k = 16;
	[V, EV, info] = eigs(W, D, k, 'lm', opts);
	EV = diag(EV);
	nof1 = sum(abs(EV-1) < 1e-6); % liczba wartości własnych bliskich "1"
	if nof1 > 1
		warning(['Wykryto ', num2str(nof1), ' wartosci wlasnych bliskich 1']);
		disp(sort(EV, 'descend'));
	end
	[m, i] = max(EV);
	y = V(:,i);
	x = D*y;
	x = x/sum(x);
otherwise % 0 < e < 1
	DW = @(x) (D*x - e*W*x);
	b = (1-e)*b;
	if exist('OCTAVE_VERSION')
	% Octave
		[y, flag, relres, iter, resid] = pcr(DW, b,  tol, maxit, D);
		fprintf(stderr, 'PCR: \n');
	else
	% MATLAB
		[y, flag, relres, iter, resid] = minres(DW, b, tol, maxit, D);
		stderr = 2; % deskryptor stderr w MATLABie
		fprintf(stderr, 'MINRES: \n');
	end
	semilogy(abs(resid));
	fprintf(stderr, '\tKod zakonczenia: %d\n\tRelRes: %e\n\tIteracji: %d\n', flag, relres, prod(iter));
	x = D*y;
end
end

Listing 7.4 - Hequation.m. [Pobierz]

%%
%% Listing 7.4 w książce "Obliczenia inżynierskie i naukowe", PWN 2011
%% (C) Piotr Krzyżanowski, 2011
%%
%% Możesz używać tego kodu zgodnie z licencją Creative Commons BY-SA, 
%% zob. http://creativecommons.org/licenses/by-sa/3.0/pl/
%%

function y = Hequation(x)
% x może być wektorem lub macierzą

global c;
global S;

y = x - 1.0./(1.0 - c*S*x);
end

Listing 7.5 - newton.m. [Pobierz]

%%
%% Listing 7.5 w książce "Obliczenia inżynierskie i naukowe", PWN 2011
%% (C) Piotr Krzyżanowski, 2011
%%
%% Możesz używać tego kodu zgodnie z licencją Creative Commons BY-SA, 
%% zob. http://creativecommons.org/licenses/by-sa/3.0/pl/
%%

function [x, resid, info, output] = newton(nazwa_f, x0, atol, rtol, step, maxit)
%
%	[x, resid, info, output] = newton(nazwa_f, x0, atol, rtol, step, maxit)
%
% Rozwiązuje układ równań nieliniowych F(x) = 0, gdzie F jest zadana M-funkcją o 
% nazwie "nazwa_f".
% Parametry wejścia:
% 	nazwa_f - nazwa M-funkcji zadającej funkcję F(). 
% 		Ta funkcja *musi* przyjmować także argumenty macierzowe!
% 	x0 - przybliżenie początkowe rozwiązania
%	atol, rtol - tolerancja błędu residualnego
%	step - co ile iteracji uaktualniać macierz pochodnej; 
%		step = 1 daje standardową metodę Newtona, 
%		step > maxit daje metodę cięciw
%	maxit - maksymalna dopuszczalna liczba iteracji
% Wyjście:
%	x - końcowe przybliżenie rozwiązania
%	resid - wektor zawierający historię zbieżnosci, 
%		resid(i) - wartość normy residuum w i-tej iteracji
%	info - kod zakończenia: = 1 - sukces, = 0 - przekroczona maxit
%	output - struktura zawierająca statystyki

% kilka pomocniczych definicji
epsq = sqrt(eps);	% pierwiastek z podwojonej precyzji maszyny
unit = ones(size(x0));	% wektor złożony z samych jedynek

% początkowe rozwiązanie przybliżone
x = x0;
	
% wyznaczenie residuum początkowego
fx = feval(nazwa_f, x);
iter = 1;
resid(iter) = norm(fx);

while( (resid(iter) >= (atol + rtol*resid(1))) & (iter < maxit) )
	if (mod(iter-1, step) == 0)
		% aproksymacja pochodnej - wszystkie kolumny jednocześnie
		h = epsq*max(abs(x),unit);
		xh = x*unit' + diag(h);
		Df = (feval(nazwa_f, xh) - fx*unit')./(unit*h');
		% rozkład macierzy pochodnej do wykorzystania
		[L, U, p] = lu(Df, 'vector');
	end
	% poprawka przybliżonego rozwiązania, x = x - Df\fx;
	x = x - U\(L\fx(p));
	
	% wyznaczenie residuum
	fx = feval(nazwa_f, x);
	iter = iter + 1;
	resid(iter) = norm(fx);
end
info = ~(iter == maxit);
output.iterations = iter - 1;
end

Listing 7.7 - solvevdw.m. [Pobierz]

%%
%% Listing 7.7 w książce "Obliczenia inżynierskie i naukowe", PWN 2011
%% (C) Piotr Krzyżanowski, 2011
%%
%% Możesz używać tego kodu zgodnie z licencją Creative Commons BY-SA, 
%% zob. http://creativecommons.org/licenses/by-sa/3.0/pl/
%%

R 	= physical_constant('MOLAR_GAS_CONSTANT'); % tylko Octave-forge; 
% w MATLAB-ie wpiszemy R = 8.314472;
a 	= 0.3640e-12; 	% stałe dla ditlenku węgla
b 	= 42.67e-6;

T 	= 303;		% temperatura
P 	= 1.0133e5;	% ciśnienie	
N 	= 0.22722;	% liczba moli gazu

V0 = N*R*T/P % objetosc gazu idealnego z rownania Clapeyrona

vdw = @(V) (P+(N^2*a)./(V.^2)).*(V-N*b)-N*R*T;

[V, FV, info, output, fjac] = fsolve(vdw, V0)

%{
% wariant z fzero
[V, FV, info, output] = fzero(vdw, V0);

% wariant z wyznaczeniem miejsc zerowych wielomianu
polyvdw = [P, -N*(b*P+R*T), N^2*a, -N^3*a*b];
V = roots(polyvdw)
%}

Listing 7.9 - vanderpoleps.m. [Pobierz]

%%
%% Listing 7.9 w książce "Obliczenia inżynierskie i naukowe", PWN 2011
%% (C) Piotr Krzyżanowski, 2011
%%
%% Możesz używać tego kodu zgodnie z licencją Creative Commons BY-SA, 
%% zob. http://creativecommons.org/licenses/by-sa/3.0/pl/
%%

function dY = vanderpoleps(t, Y, epsilon)
dY = [Y(2); ...
        epsilon*(1-Y(1)^2)*Y(2) - Y(1)];
end

Listing 7.10 - vanderpoltest.m. [Pobierz]

%%
%% Listing 7.10 w książce "Obliczenia inżynierskie i naukowe", PWN 2011
%% (C) Piotr Krzyżanowski, 2011
%%
%% Możesz używać tego kodu zgodnie z licencją Creative Commons BY-SA, 
%% zob. http://creativecommons.org/licenses/by-sa/3.0/pl/
%%

% oscylator Van der Pola

% parametry zadania
epsilon = 0.8;
Tmin = 0; 
Tmax = 60; 
Y0 = [0.01; 0]; % warunek początkowy

% t - zestaw $N$ punktów, w których chcemy wyznaczyć rozwiązanie y na przedziale [Tmin, Tmax]
N = 1000;
t = linspace (Tmin, Tmax, N);

% rozwiąż! Y - tablica wartości rozwiązania w punktach t

if exist('OCTAVE_VERSION') 
% wersja dla Octave
	vanderpol = @(Y,t) vanderpoleps(t,Y,epsilon);
	[Y, info, msg] = lsode(vanderpol, Y0, t);
	if (info ~= 2)
		disp(['Problemy w dzialaniu lsode: ', msg]);
	end
else
% wersja dla MATLAB-a
 	vanderpol = @(t,Y) vanderpoleps(t,Y,epsilon);
	[t, Y] = ode45(vanderpol, t, Y0);
end

% wizualizacja i zapis wykresów do pliku
filename = ['vanderpol', 'T-', num2str(Tmax,'%g'), ...
		'epsilon-', num2str(epsilon,'%g')]; % główny człon nazwy pliku

plot(t, Y(:,1)); % wykres rozwiązania $y(t)$
title('Oscylator Van der Pola');
xlabel('czas (t)'); 
legend('y(t)');
print('-depsc', [filename, '-sol.eps']);

input('Nacisnij dowolny klawisz...'); % czekaj na naciśnięcie klawisza!

plot(Y(:,1), Y(:,2)); % wykres fazowy rozwiązania $(y(t),y'(t))$
title('Krzywa fazowa oscylatora Van der Pola');
xlabel('y_1(t)'); ylabel('y_2(t)');
print('-depsc', [filename, '-phase.eps']);

Listing 7.11 - lorenztest.m. [Pobierz]

%%
%% Listing 7.11 w książce "Obliczenia inżynierskie i naukowe", PWN 2011
%% (C) Piotr Krzyżanowski, 2011
%%
%% Możesz używać tego kodu zgodnie z licencją Creative Commons BY-SA, 
%% zob. http://creativecommons.org/licenses/by-sa/3.0/pl/
%%

% rozwiązanie równania Lorenza

sigma = 10;
r = 28;
b = 8/3;

% t - zestaw punktów, w których chcemy wyznaczyć rozwiązanie y 
% na przedziale [0,T]

T = 200;
N = 10000; 	% wybieramy bardzo dużo punktów, by mieć lepszą wizualizację skomplikowanego rozwiązania
t = linspace (0, T, N);

% y0 - warunek początkowy

y0 = [0.01; 0; 0];

% rozwiąż! y - tablica wartości rozwiązania w punktach t

if exist('OCTAVE_VERSION') 
% wersja dla Octave
	[y, info, msg] = lsode( @(y,t)(lorenz(t,y,sigma,r,b)), y0, t);
	msg
else
% wersja dla MATLABa
	[t, y] = ode15s( @(t,y)(lorenz(t,y,sigma,r,b)), y0, t);
end
% wizualizacja i zapis do pliku

plot(y(:,1), y(:,2), 'r-');
title('Rozwiazanie modelu Lorenza: rzut 2D');
xlabel('x(t)'); ylabel('y(t)');

print('-depsc', 'lorenz2d.eps');

plot3(y(:,1), y(:,2), y(:,3), 'b-');
title('Rozwiazanie modelu Lorenza: rzut 3D');
xlabel('x(t)'); ylabel('y(t)'); zlabel('z(t)');

print('-depsc', 'lorenz3d.eps');

Listing 7.12 - lorenz.m. [Pobierz]

%%
%% Listing 7.12 w książce "Obliczenia inżynierskie i naukowe", PWN 2011
%% (C) Piotr Krzyżanowski, 2011
%%
%% Możesz używać tego kodu zgodnie z licencją Creative Commons BY-SA, 
%% zob. http://creativecommons.org/licenses/by-sa/3.0/pl/
%%

function ydot = lorenz(t,y,sigma,r,b)
ydot = [sigma*(y(2)-y(1)); ... 
	-y(1)*y(3) + r*y(1)-y(2); ...
	y(1)*y(2) - b*y(3)];
end

Listing 7.13 - cvrdrhs.m. [Pobierz]

%%
%% Listing 7.13 w książce "Obliczenia inżynierskie i naukowe", PWN 2011
%% (C) Piotr Krzyżanowski, 2011
%%
%% Możesz używać tego kodu zgodnie z licencją Creative Commons BY-SA, 
%% zob. http://creativecommons.org/licenses/by-sa/3.0/pl/
%%

function [du, flag, data] = cvrdrhs(t, u, data)
global Lap2D;
global beta;

du = -beta*(Lap2D*u) + u - u.^3;
flag = 0;
end

Listing 7.14 - lap1ddset.m. [Pobierz]

%%
%% Listing 7.14 w książce "Obliczenia inżynierskie i naukowe", PWN 2011
%% (C) Piotr Krzyżanowski, 2011
%%
%% Możesz używać tego kodu zgodnie z licencją Creative Commons BY-SA, 
%% zob. http://creativecommons.org/licenses/by-sa/3.0/pl/
%%

function L = lap1ddset(N, h)
% tworzy macierz rzadką, odpowiadajacą macierzy 
% (minus!) jednowymiarowego laplasjanu, z zerowym warunkiem Dirichleta
% N - liczba węzłów wewnętrznych, h - krok siatki

ii = [1:N, 2:N,   1:N-1];
jj = [1:N, 1:N-1, 2:N];
vv = [2*ones(1,N), -ones(1,2*(N-1))]/(h^2);
L = sparse(ii,jj,vv,N,N);
end

Listing 7.15 - cvrunrd.m. [Pobierz]

%%
%% Listing 7.15 w książce "Obliczenia inżynierskie i naukowe", PWN 2011
%% (C) Piotr Krzyżanowski, 2011
%%
%% Możesz używać tego kodu zgodnie z licencją Creative Commons BY-SA, 
%% zob. http://creativecommons.org/licenses/by-sa/3.0/pl/
%%

% Należy odkomentować linię z oznaczeniem wybranej metody 

SOLVER = 'CVODES';	% MATLAB/OCTAVE, o ile zainstalowano sundialsTB
% SOLVER = 'ODE15S';	% tylko MATLAB
% SOLVER = 'ODE5R';	% tylko Octave

if ~exist('OCTAVE_VERSION')
% MATLAB
	stderr = 2;
	stb = '/home/przykry/Matlab/sundialsTB';
else
% Octave
	stb = '/home/przykry/Octave/sundialsTB';
end

addpath(stb);
startup_STB(stb); % inicjalizacja sundialsTB

global Lap2D;
global beta;

%------------------------
% PARAMETRY i PODSTAWOWE OBIEKTY
%------------------------

% parametry zadania

a = 1; b = 1; % rozmiary obszaru $\Omega$
Tmin = 0.0; Tmax = 15.0; % długość przedziału czasowego
beta = 5e-4; % współczynnik dyfuzji

% parametry dyskretyzacji przestrzennej

Nx = 266; Ny = 267; % wewnętrzne węzły siatki
NEQ = Nx*Ny; % rozmiar zadania
hx = a/(Nx+1); hy = b/(Ny+1); % rozmiar oczka siatki
xx = hx*[1:Nx]'; yy = hy*[1:Ny]'; % współrzędne węzłów dyskretyzacji 

% dyskretny operator dyfuzji

Lap1Dx = lap1ddset(Nx, hx); Ix = speye(Nx, Nx);
Lap1Dy = lap1ddset(Ny, hy); Iy = speye(Ny, Ny);
Lap2D = kron(Iy, Lap1Dx) + kron(Lap1Dy, Ix);

clear Lap1Dx;
clear Lap1Dy;

%------------------------
% DANE POCZĄTKOWE
%------------------------

u0 = zeros(Nx,Ny);

% dane początkowe wczytujemy z pliku w formacie row-major-order (jak w OpenDX)
u0datNx = 256; u0datNy = u0datNx;
u0dat = load(['cvrd-N', num2str(u0datNx), '-step0.txt']);
u0dat = permute(u0dat,[2 1]); % indeks wiersza ma odpowiadać współrzędnej OX
mNx = min(Nx,u0datNx); mNy = min(Ny,u0datNy);
u0(1:mNx,1:mNy) = u0dat(1:mNx,1:mNy);

% musimy "rozprostować" wektor tak, by mogły na nim działać solvery ODE
u0 = u0(:);

%------------------------
% PARAMETRY PRACY SOLVERA
%------------------------

% toleracja błędu
rtol = 1.0e-5; atol = 1.0e-6;

switch SOLVER

case 'CVODES'
	data.P = []; % przestrzeń robocza
	options = CVodeSetOptions('UserData', data, 'RelTol', rtol, 'AbsTol', atol, 'LinearSolver', 'GMRES');
	CVodeInit(@cvrdrhs, 'BDF', 'Newton', Tmin, u0, options);
case {'ODE15S', 'ODE5R'}
	options = odeset('AbsTol', atol, 'RelTol', rtol, 'JPattern', (Lap2D ~= 0), 'Stats', 'on'); 
otherwise
	error('Wybrano nieznana metode');
end

%------------------------
% ROZWIĄZANIE RÓWNANIA RÓŻNICZKOWEGO
%------------------------

nout = 100;
T = linspace(Tmin, Tmax, nout+1); % chwile, w których chcemy poznać rozwiązanie

U = zeros(nout+1, NEQ); % tablica na zapisanie wartości rozwiązania w chwilach T
U(1,:) = u0'; % profil początkowy

fprintf(stderr,'Start symulacji\n'); tic;
switch SOLVER
	case 'CVODES'
		for i = 2:nout+1
			tout = T(i);
			[status,t,u] = CVode(tout,'Normal');
			U(i,:) = u'; T(i) = t;
			si = CVodeGetStats;
			fprintf(stderr, 'status = %d   t = %.2e   nst = %d   q = %d   h = %.2e\n', ...
			 status, t, si.nst, si.qlast, si.hlast);
		end
		CVodeFree;
	case 'ODE15S'
		[T, U] = ode15s(@cvrdrhs, T, u0, options); 
	case 'ODE5R'
		[T, U] = ode5r(@cvrdrhs, [Tmin, Tmax], u0, options);
		nout = length(T)-1;
end
wct = toc; % pomiar czasu działania solvera

%------------------------
% ZAPIS lub WIZUALIZACJA WYNIKÓW
%------------------------

fprintf(stderr,'Czas obliczen: %g (s)\n', wct);
filename = [lower(SOLVER),'-',num2str(Nx),'x',num2str(Ny),'-step'];
for i = 1:nout+1
% animacja rozwiązania w czasie
	t = T(i); u = reshape(U(i,:),Nx,Ny);
	imagesc(xx,yy,u);
	axis('xy', 'square');	
	title(['t = ',num2str(t)]);
	if (i == 1) || (i == nout+1) 			
		colorbar('EastOutside');
		print('-depsc',[filename,num2str(i-1),'.eps']);
	end
	pause(0.1);
end

Listing 8.1 - brus.m. [Pobierz]

%%
%% Listing 8.1 w książce "Obliczenia inżynierskie i naukowe", PWN 2011
%% (C) Piotr Krzyżanowski, 2011
%%
%% Możesz używać tego kodu zgodnie z licencją Creative Commons BY-SA, 
%% zob. http://creativecommons.org/licenses/by-sa/3.0/pl/
%%

function dx = brus(x, t)
dx = [1+x(1)^2*x(2)-4*x(1); ...
	3*x(1)-x(1)^2*x(2)];
end

Listing 8.2 - brusselator.cc. [Pobierz]

/*
**
** Listing 8.2 w książce "Obliczenia inżynierskie i naukowe", PWN 2011
** (C) Piotr Krzyżanowski, 2011
**
** Możesz używać tego kodu zgodnie z licencją Creative Commons BY-SA, 
** zob. http://creativecommons.org/licenses/by-sa/3.0/pl/
**
*/

/* Moduł Octave */
/* kompilacja: mkoctfile brusselator.cc */

#include <octave/oct.h>

DEFUN_DLD (brusselator, args, , "Funkcja prawej strony brusselatora")
{
  
  ColumnVector dx(2);
  ColumnVector x = (ColumnVector) args(0).vector_value();

  dx(0) = 1.0 + x(0)*x(0)*x(1) - 4.0*x(0);
  dx(1) = 3.0*x(0) - x(0)*x(0)*x(1);

  return octave_value (dx);
}

Listing 8.4 - brusselator.c. [Pobierz]

/*
**
** Listing 8.4 w książce "Obliczenia inżynierskie i naukowe", PWN 2011
** (C) Piotr Krzyżanowski, 2011
**
** Możesz używać tego kodu zgodnie z licencją Creative Commons BY-SA, 
** zob. http://creativecommons.org/licenses/by-sa/3.0/pl/
**
*/

/* brusselator - wersja MEX */

#include "mex.h"

void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] )
{
double *T; /* wskaźnik do skalara wejściowego */
double *X; /* wektor wejściowy */
double *dX; /* wektor wyjściowy, F(x) */

/* Check for proper number of arguments. */
if(nrhs!=2) {
	mexErrMsgTxt("Wymagane parametry: T (skalar) i X (wektor kolumnowy o 2 wspolrzednych");
} 
else if(nlhs>1) {
	mexErrMsgTxt("Funkcja zwraca tylko jeden wynik");
}

/* Create matrix for the return argument. */
plhs[0] = mxCreateDoubleMatrix(2,1, mxREAL);

/* Assign values/pointers to each input and output. */
T = mxGetPr(prhs[0]); 
X = mxGetPr(prhs[1]); 
dX = mxGetPr(plhs[0]);

dX[0] = 1.0 + X[0]*X[0]*X[1] - 4.0*X[0];
dX[1] = 3.0*X[0] - X[0]*X[0]*X[1];
}