solution_hw11 - function [f]=NewtonCoates(fun,xl,xu,n,tp)...

Info iconThis preview shows page 1. Sign up to view the full content.

View Full Document Right Arrow Icon
This is the end of the preview. Sign up to access the rest of the document.

Unformatted text preview: function [f]=NewtonCoates(fun,xl,xu,n,tp) h=(xu-xl)/n; x=xl:h:xu; for i=1:length(x) fx(i)=feval(fun,x(i)); end f=0; if (tp==1)%trapezoid integration f=fx(1)+fx(n+1); for i=2:n f=f+fx(i)*2; end f=f*h/2; elseif (tp==2) % Simpson rule if (mod(n,2)==0) %1/3 rule f=fx(1)+fx(n+1); for i=2:n if (mod(i,2)==0) f=f+fx(i)*4; else f=f+fx(i)*2; end end f=f*h/3; else % combined 1/3 and 3/8 rule m=n-3; f=fx(1)+fx(m+1); for i=2:m if (mod(i,2)==0) f=f+fx(i)*4; else f=f+fx(i)*2; end end f=f*h/3; % the 1/3 rule part end f=f+3*h/8*(fx(n+1)+3*fx(n)+3*fx(n-1)+fx(n-2)); % the 3/8 part end -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.function [f]=test(x) f=2/(1+2*x^2); We can also write a program to carry out the comparison required in the homework, but if you did it manually and didn’t present the figure in your homework, there is also no penalty. segnum=[4 5 6 7 100 1000]; xl=-3; xu=3; for i=1:length(segnum) n=segnum(i); int_Trap(i)=NewtonCoates('test',xl,xu,n,1); int_Simpson(i)=NewtonCoates('test',xl,xu,n,2); end figure(1) semilogx(segnum,int_Newton,'o-',segnum,int_Simpson,'*-' ); grid on; legend('Trapzoid rule integration','Simpson rule intergration') xlabel('segment number'); ylabel('integration'); ylim([2 6]) 6 5.5 5 4.5 integration 4 3.5 3 2.5 2 0 10 Trapzoid rule integration Simpson rule intergration 10 10 segment number 1 2 10 3 Comparison table: Number of Integration by Integration by segment Trapezoidal rule Simpson rule 4.2488 3.5598 4 3.5587 3.4364 5 3.8830 4.2183 6 3.7306 3.7507 7 3.7896 3.8666 10 3.7881 3.7882 100 3.7882 3.7882 1000 ...
View Full Document

Ask a homework question - tutors are online