《计算方法》上机实验试题
注:以下代码均通过matlab语言实现
1. (25分)计算积分
若用下列两种算法
试依据积分In的性质及数值结果说明何种算法更合理。
测试代码question1.m:
format short e
% 算法A:由1-20开始递推
I_forward=zeros(21,1);
I_forward(1)=log(6/5); % 先计算I0,但MATLAB中下标必须从1开始
for n=2:21
I_forward(n)=-5*I_forward(n-1)+1/(n-1); % 课本上1/n要改为1/(n-1)
end
fprintf('In calculated from n=1 to n=20:\n');
I_forward(2:21)
% 算法B:由20-1开始递推
I_backward=zeros(21,1);
I_backward(21)=0.5*(1/(6*21)+1/(5*21)); % 由放缩不等式计算I20,这里下标为21
for n=21:-1:2
I_backward(n-1)=-(1/5)*I_backward(n)+1/(5*(n-1)); % 1/n同样要改为1/(n-1)
end
fprintf('In calculated from n=19 to n=0:\n');
I_backward(20:-1:1)
计算结果(这里每道题结果均为MATLAB截图):
左图对应算法A,从上到下为I1-I20(I0为手工计算);右图对应算法B,从上到下为I19-I0(I20为手工计算)。
性质:In>0且单调递减。算法A从n=18开始并未出现单调递减,其每一步舍入误差都增长5倍;而算法B每一步舍入误差都减少1/5。所以算法B合理。
2. (25分)求解方程f(x)=0有如下牛顿迭代公式
(1) 编制上述算法的通用程序,并以 (ε为预定精度)作为终止迭代的准则。
(2) 利用牛顿迭代求解方程,设预定精度ε=10-10。
向量型牛顿迭代法的通用程序Newton_iteration.m,传入的参数均可以是向量:
% newton_iteration: 向量型Newton迭代法
function [x1,k] = newton_iteration(f,f_de,varible,x0,eplison)
% dx 表示微分,f_de表示f的导数,variable表示f中的变量
% subs函数将f和f_de中的变量替换为x0中的值
dx = -subs(f_de,varible,x0)\subs(f,varible,x0);% A\b相当于inv(A)*B
x1=x0+dx;
k=1; % 迭代次数
while norm(x1-x0)>=eplison
x0=x1;
dx=-subs(f_de,varible,x0)\subs(f,varible,x0);
x1=x0+dx;
k=k+1;
end
end
测试代码question2.m:
format short
syms x y
f1=5*cos(x)-x*y-3;
f2=y^2+8*x*y-7;
F=[f1;f2];
F_de=[diff(F,'x'),diff(F,'y')]; % 这里生成4*4的雅克比矩阵
eplison=10e-10;
variable=[x;y];
X0=[1;1];
[x,k]=newton_iteration(F,F_de,variable,X0,eplison);
fprintf('after %d newton_iteration, the root is:\n', k);
for i=1:length(x)
fprintf('%1.8f\n',x(i));
end
计算结果:
由于eplison=10e-10计算时间过长,上面是将精度改成10e-4时的计算结果。
3. (25分) 已知,
(1) 利用插值节点x0=1.00,x1=1.02,x2=1.04,x3=1.06,构造三次Lagrange插值公式,由此计算f(1.03)的近似值,并给出其实际误差;
(2) 利用插值节点x0=1,x1=1.05构造三次Hermite插值公式,由此计算f(1.03)的近似值,并给出其实际误差。
使用拉格朗日插值算法的函数lagrange.m:
%% lagrange: Lagrange插值函数
function Y = lagrange(X_inter,Y_inter,X)
% X_inter是插值节点,Y_inter是插值节点的函数值,X是待求点,Y是待求点的函数值
n=length(X_inter); % 插值节点个数
m=length(X); % 待求节点个数
% 下面这个循环求每个待求点的函数值
for X_num=1:m
x=X(X_num);
P=0.0; % P是Lagrange插值函数
% 下面这个循环求Lagrange插值函数P,即所有基函数L_j(x)*y_j的和
for j=1:n
L=1.0; % L是Lagrange基函数
% 下面这个循环求Lagrange基函数L_j(x)
for i=1:n
if i~=j
L=L*(x-X_inter(i))/(X_inter(j)-X_inter(i));
end
end
P=P+L*Y_inter(j);
end
Y(X_num)=P;
end
end
使用埃尔米特插值算法的函数hermite_interpolation.m:
%% hermite_interpolation: Hermite插值算法
function Y = hermite_interpolation(X_inter,Y_inter,Y_inter_de,X)
% Y_inter_de代表插值节点的导数值
n=length(X_inter); % 插值节点个数
m=length(X); % 待求节点个数
% 下面这个循环求每个待求点的函数值
for X_num=1:m
x=X(X_num);
H=0.0; % H是hermite插值函数
% 下面这个循环求hermite插值函数H
for j=1:n
L_2=1.0; % L_2是Lagrange基函数的平方
b=0; % b是求解H的中间量
% 下面这个循环求Lagrange基函数L_j(x)的平方
for i=1:n
if i~=j
L_2=L_2*((x-X_inter(i))/(X_inter(j)-X_inter(i)))^2;
b=1/(X_inter(j)-X_inter(i))+b;
end
end
H=H+L_2*((X_inter(j)-x)*(2*b*Y_inter(j)-Y_inter_de(j))+Y_inter(j));
end
Y(X_num)=H;
end
end
测试代码question3.m:
% 第一小题,利用Lagrange插值公式
syms x
f(x)=exp(x)*(3*x-exp(x));
X_inter_lag=[1.00;1.02;1.04;1.06];
Y_inter_lag=f(X_inter_lag);
X=1.03;
Y_lag=lagrange(X_inter_lag,Y_inter_lag,X);
error_lag=abs(Y_lag-f(X)); % 实际误差
fprintf('P(1.03)=%1.8f(cal by lagrange)\n', Y_lag);
fprintf('error=%1.8f(cal by lagrange)\n', error_lag);
% 第二小题,利用Hermite插值公式
X_inter_her=[1;1.05];
Y_inter_her=f(X_inter_her);
f_de=diff(f);
Y_inter_her_de=f_de(X_inter_her);
Y_her=hermite_interpolation(X_inter_her,Y_inter_her,Y_inter_her_de,X);
error_her=abs(Y_her-f(X)); % 实际误差
fprintf('H(1.03)=%1.8f(cal by hermite_interpolation)\n', Y_her);
fprintf('error=%1.8f(cal by hermite_interpolation)\n', error_her);
计算结果:
(25分) 利用复合求积算法计算积分,取步长为10-4
使用复合梯形公式的函数compound_trapezoidal.m:
%% compound_trapezoidal: 复合梯形公式
function F = compound_trapezoidal(f,a,b,h)
% f是被积函数,a,b是积分区间,h是步长
N=(b-a)/h; % 节点个数
f_sum=0;
for n=1:N-1
x_n=a+n*h;
f_sum=f_sum+f(x_n);
end
F=h/2*(f(a)+2*f_sum+f(b));
end
使用复合辛普森公式的函数compound_simpson.m:
%% compound_simpson: 复合辛普森公式
function F = compound_simpson(f,a,b,h)
% f是被积函数,a,b是积分区间,h是步长
N=(b-a)/h; % 节点个数
f_sum_half=0;
f_sum=0;
x_n=a; % 保证x_n_half的起始计算
for n=1:N-1
x_n_half=x_n+h/2;
x_n=a+n*h;
f_sum_half=f_sum_half+f(x_n_half);
f_sum=f_sum+f(x_n);
end
F=h/6*(f(a)+4*f_sum_half+2*f_sum+f(b)); % 注意课本上x_n_half的求和是从n=0开始的
end
测试代码question4.m:
syms x
g(x)=sqrt(1+cos(x)^2);
a=0;
b=48;
% h=10e-4;
h=0.1;
% 使用复合梯形公式
F_compound_trapezoidal=compound_trapezoidal(g,a,b,h);
fprintf('cal by compound_trapezoidal: %1.8f\n', F_compound_trapezoidal);
int_res=int(g,[0,48]);
fprintf('cal by int: %1.8f\n', int_res); % 使用matlab自带的求积分函数
% 使用复合辛普森公式
F_compound_simpson=compound_simpson(g,a,b,h);
fprintf('cal by compound_simpson: %1.8f\n', F_compound_simpson);
计算结果:
其中cal by int 是使用matlab自带的求积分函数int计算出的结果。