数值计算方法第3、4次实验--龙贝格--龙格库塔

1.完成复合梯形、复合辛普森求积公式，用变步长、事后误差估计的方法完成P145例题6.3.1对应表6-3

的计算，课后作业题7.

%复合梯形求积公式的MATLAB实现%姓名：王定

%学号：1306034248

%中北大学仪器与电子学院

function T_n=ComTrapezoidal(a,b,n)

%a为积分上限

%b为积分下限

%n为划分区间的个数

h=(b-a)/n;

for(k=0:n)

x(k+1)=a+k*h;

if(x(k+1)==0)

x(k+1)=10^(-10);

end

end

I_1=h/2*(f(x(1))+f(x(n+1)));

for(i=2:n)

F(i)=h*f(x(i));

end

I_2=sum(F);

I_n=I_1+I_2

%****************************% %复合辛普森求积公式的MATLAB实现

function I=ComSimpson(a,b,n)

n=2;

h=(b-a)/2;

I1=0;

I2=(f(a)+f(b))/h;

eps=1.0e-4;

while(abs(I2-I1)>eps)

n=n+1;

h=(b-a)/n;

I1=I2;

I2=0;

for(i=0:n-1)

x=a+h*i;

x1=x+h;

I2=I2+h/6*(f(x)+4*f((x+x1)/2)+f(x1));

end

end

I=I2 %变步长梯形求积法的MATLAB实

现

function AdaptiveTrapezoidal(a,b,eps)

m=1

t0=0;

t2=(f(a)+f(b))/4+f((a+b)/2)

while(abs(t2-t0)>eps)

m=m+1

p=t2;

t1=0;

n=2^m;

h=(b-a)/n;

for(k=0:(n-1))

t1=t1+h*((f(a+k*h)+f(a+(k+1)*h))/4+

f(a+(k+1/2)*h)/2);

end

t2=t1

t0=p;

end

I=t2

%**************************%

>>

AdaptiveTrapezoidal(1,2,0.000001)

m = 1

t2 = 3.03948481584447

m = 2

t2 = 2.02304986763725

m = 3

t2 = 2.02080858246806

m = 4

t2 = 2.02024624995310

m = 5

t2 = 2.02010553934348

m = 6

t2 = 2.02007035369492

m = 7

t2 = 2.02006155678257

m = 8

t2 = 2.02005935752321

m = 9

t2 = 2.02005880770641

I = 2.02005880770641

%课后作业题7.

>>

AdaptiveTrapezoidal(1,

3,0.0001)

m = 1

t2 =

7.99930630234471

m = 2

t2 =

10.84204346755743

m = 3

t2 =

10.92309388961378

m = 4

t2 =

10.94339842118680

m = 5

t2 =

10.94847716722413

m = 6

t2 =

10.94974701694155

m = 7

t2 =

10.95006448956964

m = 8

t2 =

10.95014385836406

I =

10.95014385836406

2.完成龙贝格积分，计算P150例题6.4.2对应的表6-5，按例题列表格显示数值，完成课后作业题8.

%龙贝格求积法的MATLAB实现

%姓名：王定

%学号：1306034248

%中北大学仪器与电子学院

function [I,step]=Romberg(f,a,b,eps)

f=input('please input f=');

a=input('please input a=');

b=input('please input b=');

eps=input('please input eps=');

%I为所求积分值

%step为积分划分的子区间次数

%f为被积函数

%a为积分上限

%b为积分下限

%eps为积分精度

if(nargin==3)

eps=1.0e-4;

end;

M=1;

tol=10;

k=0;

T=zeros(1,1);

h=b-a;

T(1,1)=(h/2)*(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsy m(sym(f)),b));

while(tol>eps)

k=k+1;

h=h/2;

Q=0;

for(i=1:M)

x=a+h*(2*i-1);

Q=Q+subs(sym(f),findsym(sym(f)),x);

end

T(k+1,1)=T(k,1)/2+h*Q;

M=2*M;

for(j=1:k)

T(k+1,j+1)=T(k+1,j)+(T(k+1,j)-T(k,j))/(4^j-1);

end

tol=abs(T(k+1,j+1)-T(k,j));

end

I=T(k+1,k+1)

step=k %计算P150例题6.4.2对应的表6-5

>> [I,step]=Romberg('4/(1+x^2)',0,1); format long

I,step

I =

3.14159266527772

step =

4

%********************************%

%按例题列表格显示数值，完成课后作业题8.

>>

please input f='(2/sqrt(pi))*exp(-x)'

please input a=0

please input b=1

please input eps=1.0e-5

I =

0.71327166981418

step =

3