Prev: I have a PROBLEM
Next: actxcontrol
From: Koubar os on 13 Feb 2006 11:09
I want to 'translate' some programs I had in C for matlab but by
sigtly optimizing the code for matlab use. I am new in matlab. So, I
had the following code:

(1) for K = 1:dT:M
..
(2) for I = 1:N+1
dU(I) = (U(I+1)-2*U(I)+U(I-1))/dX^2;
end
..
(3) for I = 1:N+1
U(I) = U(I)+dU(I)*dT
end
..
end

This actually is a simple implementation of an explicit scheme for
the numerical solution of the heat equation. Now I managed to do this
and eliminate the (2) and (3) loop:

(1) for K = 1:dT:M
..
(2) dU(1:N+1) = (U([2:end 1])-2*U+U([end 1:end-1]))/dX^2;
..
U(I+1)-->U([2:end 1])
U(I-1)-->U([end 1:end-1])
..
(3) U(1:N+1) = U(1:N+1)+dU(1:N+1)*dT
..
end

I could really use your help in order to eliminate the (1) time loop,
in a way that the values of U(1:N+1) are stored for each time step dT
or, in case M is very large this would take a lot of memory and time,
get 'snapshots' every a couple of timesteps which I could define
explicitly.
From: Greg von Winckel on 14 Feb 2006 22:49
I don't know if this will help you, but here is some old code I have
lying around from a homework problem:

function [x,u]=heateq(M,N,T,r)

% Solve the Heat Equation using 2nd order finite differences for
% spatial discretization and Crank-Nicholson for time stepping

M1=M+1; % # of grid points = # of intervals + 1
x=linspace(-1,1,M1)'; % Set up grid points
h=x(2)-x(1); % Grid spacing
k=T/N; % Time step size
del=k/(2*h^2); % Crank-Nicholson parameter
in=2:M; m=M-1; % indices of interior points
u=sin(pi*r*x); % Get initial condition
e=ones(M-1,1); % Column vector of all ones

% Left hand side matrix
A=spdiags([-del*e, (1+2*del)*e, -del*e], -1:1, m, m);

[L,U]=trilu(A); % Do LU factorization
coeff=1-2*del; % Lump diagonal coefficient for RHS

for n=1:N % Time stepping loop

f=coeff*u+del*([0;u(1:M)]+[u(2:M1);0]); %RHS
u(in)=U\(L\f(in));
end

exact=exp(-(pi*r)^2*T)*sin(pi*r*x);
err=norm(u-exact,1)/norm(exact,1);

function [L,U] = trilu(A);
n=length(A); am=diag(A,-1); a0=diag(A); ap=diag(A,1);
for k = 1:(n-1)
am(k)=am(k)/a0(k); k1=k+1; a0(k1)=a0(k1)-am(k)*ap(k);
end
L=speye(n)+spdiags(am,-1,n,n); U=spdiags([a0 [0;ap]],0:1,n,n);
 | 
Pages: 1
Prev: I have a PROBLEM
Next: actxcontrol