hw10_p03 - end end e%Make the F vector%All F's need...

Info iconThis preview shows pages 1–2. Sign up to view the full content.

View Full Document Right Arrow Icon
clc clear D = input('D: '); K = input('K: '); h = input('h: '); xmin = input('x_min: '); xmax = input('x_max: '); Y(1) = input('y at x_min: '); x = xmin:h:xmax; n = length(x); Y(n) = input('y at x_max: '); int_1 = 1; i %For this linear BVP, we have: %p(x) = 0 (no coefficient for y') %q(x) = K/D %f(x) = 0 % %We need to solve for the y values that aren't endpoints (known values), %meaning length(x)-2 y's. So our matrix will have dimensions length(x)-2 %by length(x)-2 % %Start with a T matrix of appropriate dimensions which is all zeros T = zeros( n-2, n-2 ); T for i = 1:n-2 p(i) = 0.0; q(i) = K/D; f(i) = 0.0; % q(i) = 1 - x(i+1); % f(i) = x(i+1); end e for i=1:n-2 f %Avoid having an array index of 0 by using an if loop if ( i > 1 ) T(i,i-1) = -(1.0 + 0.5*p(i)*h); end T(i,i) = 2.0 + q(i)*h^2; %Avoid having an array index > length(x)-2 by using an if loop if ( i < n-2 ) T(i,i+1) = -(1.0 + 0.5*p(i)*h);
Background image of page 1

Info iconThis preview has intentionally blurred sections. Sign up to view the full version.

View Full DocumentRight Arrow Icon
Background image of page 2
This is the end of the preview. Sign up to access the rest of the document.

Unformatted text preview: end end e %Make the F vector: % %All F's need h^2*f(xi) term: for i=1:n-2 if (i == 1) F(i,1) = (1.0 + 0.5*p(i)*h)*Y(1) + (h^2)*f(i); elseif (i == n-2) F(i,1) = (1.0 - 0.5*p(i)*h)*Y(n) + (h^2)*f(i); else F(i,1) = (h^2)*f(i); end end e U = T\F; U for i=1:n-2 Y(i+1) = U(i,1); end e % open file fid4 = fopen('hw10_p03_answer.txt','a'); % 'wt' means &amp;quot;write text&amp;quot; if (fid4 &amp;lt; 0) error('could not open file &amp;quot;hw10_p03_answer.txt&amp;quot;'); end; fprintf(fid4,'\n\nh = %10.5f\n', h); fprintf(fid4,'y(%10.5f)=%10.5f, y(%10.5f)=%10.5f\n', xmin, Y(1), xmax, Y(n)); fprintf(fid4,'\t\t\tx\t\t\t\tY1(x)\n'); for i = 1:n if mod(i-1,int_1) == 0 fprintf(fid4,'%15.5f%15.5f\n', x(i),Y(i)); end end fclose(fid4); f figure(1), plot(x,Y,'o',x,Y) xlabel('x (cm)'); ylabel('Concentraion (M)'); legend('C_A'); title('Problem 03 in HW#10: figure');...
View Full Document

Page1 / 2

hw10_p03 - end end e%Make the F vector%All F's need...

This preview shows document pages 1 - 2. Sign up to view the full document.

View Full Document Right Arrow Icon
Ask a homework question - tutors are online