function [T,timeV,spaceV]=nma_PDE_parabolic_explicit_rod_with_rate_BC(length,leftTRate,rightTRate,rodT,deltaT,nSteps,deltaX,k) %solve parabolic PDE using explicit method for changing boundary conditions %function [T,timeV,spaceV]=nma_PDE_parabolic_explicit_rod_with_rate_BC(length,leftTRate,rightTRate,rodT,deltaT,deltaX,k) % %Solves the parabolic PDE using explicit method for given Boundary %expressed as rate (dT/dt). i.e. B.C. is not fixed values. % % %INPUT: % length: length of rod in some length units % leftTRate: rate of change of independent variable (V or T for example) at % left end of rod w.r.t time % rightTRate: rate of change of independent variable (V or T for example) at % right end of rod w.r.t. time % rodT: Value of independent variable (V or T for example) at % rod interior points at t=0 (initial) % deltaT: time step size % nSteps: number of time steps to run for. % deltaX: step size in distance along rod % k: material coefficent. depends on problem (heat, Voltage, etc...) % % OUTPUT: % T : a matrix that represent the solution. Each row is the solution % at different time. i.e. first row is the solution at t=0; % second row is the solution after one time step, etc... % first column is solution at left edge, second column is solution % at one distance step to the right of the left edge, etc... % until last column is reached, which is solution at right edge. % timeV : a vector that has the time values corresponding to when each % row in T was generated. % % spaceV : a vector that has the space values corresponding to the points % on the rod. These are the points where each column in T is % located. %Author: Nasser Abbasi nIntervals=length/deltaX; nPoints=nIntervals+1; T=zeros(nSteps,nPoints); %set initial conditions; T(1,:)=rodT; timeV=zeros(nSteps,1); spaceV=zeros(nPoints,1); for(i=1:nPoints) spaceV(i)=(i-1)*deltaX; end lambda=k*deltaT/(deltaX)^2; for(t=1:nSteps-1) for(i=2:nPoints-1) %do left edge T(t+1,1)=T(t,1)+leftTRate*deltaT; T(t+1,end)=T(t,end)+rightTRate*deltaT; T(t+1,i)=T(t,i)+lambda*(T(t,i+1)-2*T(t,i)+T(t,i-1)); timeV(t+1)=timeV(t)+deltaT; end end