复化中矩形，复化梯形，复化Simpson，龙贝格函数，quad函数，integral函数
第一题
1）复化中矩形
%复化中矩形
clc,clear;
format long;
%a = 1;
%b = 2;
a = 2;
b = 3;
n = 11000;
%n = 115300;
d = (b-a)/n;
sum = 0;
%y = @(x)2*x^3*exp(x^2)/3;
y = @(x)2*x/(x^2-3);
for i = a:d:b-d
sum = sum + d*y((i+i+d)/2);
end
%disp(exp(4));
disp(sum);
disp(log(6));
%exp(4)-sum
log(6)-sum
2）复化梯形
%复化梯形
clc, clear;
a = 1;
b = 2;
%a = 2;
%b = 3;
n = 165000;
%n = 15000;
d = (b-a)/n;
sum = 0;
y = @(x)2*x^3*exp(x^2)/3;
%y = @(x)2*x/(x^2-3);
for i = a:d:b-d
sum = sum + d * (y(i)+y(i+d))/2;
end
disp(sum);
disp(exp(4));
%disp(log(6));
sum-exp(4)
%sum-log(6)
3)复化Simpson
%复化Simpson
clc, clear;
%a = 1;
%b = 2;
a = 2;
b = 3;
n = 95;
%n = 249;
d = (b-a)/n;
sum = 0;
%y = @(x)2*x^3*exp(x^2)/3;
y = @(x)2*x/(x^2-3);
for i = a:d:b-d
x0 = i;
x1 = (i + i + d)/2;
x2 = i + d;
sum = sum + d * (y(x0)+4*y(x1)+y(x2))/6;
end
disp(sum);
%disp(exp(4));
disp(log(6));
%sum-exp(4)
sum-log(6)
4)龙贝格
%龙贝格
clc, clear;
%a = 1;
%b = 2;
a = 2;
b = 3;
pre = 0.5e-8;
T = zeros(6);
%y = @(x)2.*x.^3.*exp(x.^2)./3;
y = @(x)2.*x./(x.^2-3);
h = @(m)(b-a)/(2^m);
T(1, 1) = h(1)*(y(a) + y(b));
for i = 2:10
for j = 1:i
if(j == 1)
T(i, 1) = T(i-1, 1)/2;
k = 0:2^(i-2)-1;
xi = a + k * h(i-2)+h(i-1);
T(i, 1) = T(i, 1) + sum(h(i-1)*y(xi));
else
T(i, j) = 4^(j-1)/(4^(j-1)-1)* T(i, j-1)-1/(4^(j-1)-1)*T(i-1, j-1);
end
end
end
%disp(exp(4));
disp(log(6));
之前一直把握不好n怎么求的，现在在搭档大佬的帮助下看到了推导的过程(ﾉ･ω･)ﾉ开心！撒花★°:.☆(￣▽￣)/$:.°★* 。
贴图走一波~
1.jpg
2.jpg
3.jpg
4.jpg
5.jpg
6.jpg
7.jpg
第二题
求I(4), I(8), I(32)
（1）龙贝格函数
clc;
clear;
syms x;
k1 = 4;
%k1 = 4;
%k1 = 4;
n = 8;
a = 0; %积分下限
b = k1 * pi; %积分上限
R = 0.5 * 10 ^(-6); %精度
T = zeros(n, n);
y = abs(sin(x)); %函数
p = subs(y, x, a);
q = subs(y, x, b);
if a == 0
p=1;
end
T(1, 1) = (b-a)/2*(p+q);
for i = 2:n
f = 0;
for j = 1:2^(i-2)
t = a+((2*j-1)/2^(i-1))*(b-a);
z = subs(y, x, t);
f = f+z;
end
T(i, 1) =0.5*T(i-1, 1)+f*(b-a)/2^(i-1);
end
for j = 2:n
n = n - 1;
for i = 1:n
T(i, j) = (4^(j-1) * T(i+1, j-1)-T(i, j-1))/(4^(j-1)-1);
end
end
T = vpa(T, 7);
T = eval(T);
n = 8;
w = ones(1, 7);
for j = 2:n
if abs(T(1, j) - T(1, j-1)) <= R
break;
end
end
a = T(1, j);
a = vpa(a, 7);
T
a
% k1 = 4
T =
1 至 5 列
6.283185000000000 2.094395000000000 0.977384400000000 9.558864000000000 8.214549999999999
3.141593000000000 1.047198000000000 9.424778000000000 8.219801000000000 8.119866000000000
1.570796000000000 8.901179000000001 8.238629000000000 8.120255999999999 8.059881000000001
7.068583000000000 8.280037999999999 8.122106000000000 8.060117000000000 8.029941000000001
7.977175000000000 8.131976000000000 8.061086000000000 8.030058000000000 0
8.093275999999999 8.065516000000001 8.030543000000000 0 0
8.072456000000001 8.032729000000000 0 0 0
8.042661000000001 0 0 0 0
6 至 8 列
8.119773000000000 8.059808000000000 8.029902000000000
8.059822000000000 8.029904000000000 0
8.029911000000000 0 0
0 0 0
0 0 0
0 0 0
0 0 0
0 0 0
a =
8.029902
% k1 = 8
T =
1 至 5 列
12.566369999999999 4.188790000000000 1.954769000000000 0.961870300000000 18.706109999999999
6.283185000000000 2.094395000000000 0.977384400000000 18.636790000000001 16.189579999999999
3.141593000000000 1.047198000000000 18.360859999999999 16.199130000000000 16.119969999999999
1.570796000000000 17.278759999999998 16.232910000000000 16.120280000000001 16.059880000000000
13.351770000000000 16.298279999999998 16.122039999999998 16.060120000000001 0
15.561650000000000 16.133050000000001 16.061080000000000 0 0
15.990200000000000 16.065580000000001 0 0 0
16.046740000000000 0 0 0 0
6 至 8 列
16.187120000000000 16.119879999999998 16.059799999999999
16.119900000000001 16.059809999999999 0
16.059819999999998 0 0
0 0 0
0 0 0
0 0 0
0 0 0
0 0 0
a =
16.0598
% k1 = 32
T =
1 至 5 列
50.265479999999997 16.755160000000000 7.819075000000000 3.847481000000000 1.916197000000000
25.132739999999998 8.377580000000000 3.909538000000000 1.923741000000000 0.958098300000000
12.566369999999999 4.188790000000000 1.954769000000000 0.961870300000000 73.387270000000001
6.283185000000000 2.094395000000000 0.977384400000000 73.104360000000000 64.039730000000006
3.141593000000000 1.047198000000000 71.977379999999997 64.075140000000005 0
1.570796000000000 67.544240000000002 64.198610000000002 0 0
51.050879999999999 64.407709999999994 0 0 0
61.068500000000000 0 0 0 0
6 至 8 列
0.957161700000000 73.475780000000000 64.027709999999999
73.458070000000006 64.028289999999998 0
64.030590000000004 0 0
0 0 0
0 0 0
0 0 0
0 0 0
0 0 0
a =
64.02771
（2）quad函数
% quad函数
clc;
k1 = 4;
k2 = 8;
k3 = 32;
f1 = quad(@(x) abs(sin(x)), 0, k1*pi);
f2 = quad(@(x) abs(sin(x)), 0, k2*pi);
f3 = quad(@(x) abs(sin(x)), 0, k3*pi);
f1
f2
f3
f1 =
8.0000
f2 =
16.0000
f3 =
64.0000
（3）integral函数
clc;
clear;
k1 = 4;
k2 = 8;
k3 = 32;
f = @(x) abs(sin(x));
a1 = integral(f, 0, k1*pi);
a2 = integral(f, 0, k2*pi);
a3 = integral(f, 0, k3*pi);
a1
a2
a3
a1 =
7.999999399161427
a2 =
16.000000563470667
a3 =
64.000001399531186
第三题
clf
h=0.05;
t=0:h:1;
x=t.^3; %f(x) = x^3;
dxdt_diff=diff(x)/h;
dxdt_grad=gradient(x)/h;
subplot(1,2,1) %第一子图
plot(t,x,'b')
hold on
plot(t,dxdt_grad,'m','LineWidth',8)
plot(t(1:end-1),dxdt_diff,'.k','MarkerSize',8)
axis([0,1,-0.1,1.1])
title('[0, 1]')
legend('x(t)','dxdt_{grad}','dxdt_{diff}','Location','North')
xlabel('t')
box off
hold off
subplot(1,2,2) %第二子图
kk=(length(t)-10):length(t);
hold on
plot(t(kk),dxdt_grad(kk),'om','MarkerSize',8)
plot(t(kk-1),dxdt_diff(kk-1),'.k','MarkerSize',8)
title('[end-10, end]')
xlabel('t'),box off
legend('dxdt_{grad}','dxdt_{diff}','Location','SouthEast')
hold off
第三题.jpg
网友评论