function [U,xv] = heat_pwlin(z,u0,tv,xv) % u = heat_pwlin(z,u0,tv) % u = heat_pwlin(z,u0,tv,xv) % find the solution of heat equation u_t = u_{xx} where initial value u0 is % continuous, piecewise linear, constant for x<=z(1) and x>=z(end) % inputs: % z: nodes for initial data (increasing) % u0: values of u0 at nodes % tv: t-values for which we want to find solution % xv: (optional) x-values for which we want to find solution % outputs: % u: array with values of solution at points xv, tv % each column corresponds to one tv value % Examples: % heat_pwlin([0 1 2 3],[1 -2 2 -1],0:.02:1); % heat_pwlin([0 1 2 3],[1 -2 2 -1],0.3,0.7) % find u for t=.3,x=.7 % x=-1:.1:1; heat_pwlin(x,1-x.^2,0:.02:1); % u0(x)=1-x^2 for x in [-1,1], 0 otherwise % needs colorcurves.m tv(tv==0)=1e-15; % use small number instead of zero to avoid division by zero if nargin==3 % if no xv provided, make reasonable choice for xv: d = z(end)-z(1); h = min(diff(z)); % h is smallest interval length xv = linspace(z(1)-d,z(end)+d,1+3*10*d/h); % spacing should be <= h/10 end [X,T] = ndgrid(xv,tv); G = @(x) x.*(erf(x/2)+1)/2 + exp(-x.^2/4)/pi^.5; % antiderivative of g(x) U = u0(1)*ones(size(X)); % U = u0(1) c = diff(u0)./diff(z); % slopes p = diff([0;c(:);0]); % jumps of slopes for j=1:length(z) U = U + p(j)*T.^.5.*G(T.^-.5.*(X-z(j))); % add (jump of slope at z_j)*V(x-z_j,t) end colorcurves(X,T,U) % g(x) = erf(x/2)/2 + 1/2 % antiderivative is G(x) = x*(erf(x/2)+1)/2 + exp(-x^2/4)/pi^.5 % U(x,t) = g(t^-.5*x) solution of heat eq for u0=Heaviside fct % S(x,t) fundamental solution, partial deriv. of U wrt x % V(x,t) = t^.5*G(t^-.5*x) is antiderivative wrt x of U % (we omit t argument of S,U,V below) % solution formula for integral over [a,b] where u0 is linear, and u0'=c=constant: % int_[a,b] S(x-y)*u0(y) dy % = int_[a,b] U(x-y)*c dy - U(x-b)*u0(b) + U(x-a)*u0(a) % = (-V(x-b)+V(x-a))*c - U(x-b)*u0(b) + U(x-a)*u0(a) % Because of continuity of u0(x) all terms with U(x-a) and U(x-b) will cancel % when we add contributions from all intervals together.