hw09_p02 - end T(i,i) = 2.0 + q(i)*h^2; %Avoid having an...

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

View Full Document Right Arrow Icon
clc clear 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 is for making table using 'mod' operator if h == pi/8 int_1 = 1; elseif h == pi/16 int_1 = 2; elseif h == pi/32 int_1 = 4; elseif h == pi/64 int_1 = 8; elseif h == pi/128 int_1 = 16; else int_1 = 1; end e %For this linear BVP, we have: %p(x) = 0 (no coefficient for y') %q(x) = 5 cos(3x) %f(x) = x^2 % %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) = 5.0*cos(3.0*x(i+1)); f(i) = x(i+1)^2; % 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);
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 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); 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('hw09_p02_answer.txt','a'); % 'wt' means "write text" if (fid4 < 0) error('could not open file "hw09_p02_answer.txt"'); 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);...
View Full Document

Page1 / 2

hw09_p02 - end T(i,i) = 2.0 + q(i)*h^2; %Avoid having an...

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