题目一
syms x;
fx=inline('54*x^6+45*x^5-102*x^4-69*x^3+35*x^2+16*x-4');
% 首先利用diff(54*x^6+45*x^5-102*x^4-69*x^3+35*x^2+16*x-4)求解函数的导数
dfx=inline('324*x^5 + 225*x^4 - 408*x^3 - 207*x^2 + 70*x + 16');
%% 绘图
x=-1.5:0.0001:1.2;
y=54.*x.^6+45.*x.^5-102.*x.^4-69.*x.^3+35.*x.^2+16.*x-4;
y2=0.*(1<=x);% 找到零点的大致范围，以便在牛顿迭代法中设置初值
plot(x,y,x,y2)
%% 调用函数分别求出五个根
[re1,step1]=NewtonDiedai(fx,dfx,-1.39);
[re2,step2]=NewtonDiedai(fx,dfx,-0.68);
[re3,step3]=NewtonDiedai(fx,dfx,0.23);
[re4,step4]=NewtonDiedai(fx,dfx,0.51);
[re5,step5]=NewtonDiedai(fx,dfx,1.18);
function [result,step] = NewtonDiedai(fx,dfx,x0)
%此方法为牛顿迭代法求解非线性方程的根
% fx为原函数 dfx为导函数 x0为初值
% result为结果，step为迭代步数
x1=x0;
step=1;
while 1
x0=x1;
x1=x0-feval(fx,x0)/feval(dfx,x0); %feval()函数执行指定的函数
if abs(x0-x1)<0.0001 %判断迭代结果精度是否达到要求
break;
end
if(step>50) %当步数大于50时认为迭代不收敛，退出
break;
end
step=step+1; %记录步数
end
result=x0;
end
题目二
%% 求导
syms x;
fx=(1-3/(4*x))^(1/3);
diff(fx) %求出函数的导数
%% 绘图
x=-3:0.0001:3;
y=(1-3./(4.*x)).^(1/3);
y2=0.*(1<=x);% 找到零点的大致范围，以便在牛顿迭代法中设置初值
plot(x,y,x,y2)
%% 调用函数求出根
fx_line=inline('(1-3/(4*x))^(1/3)');
dfx_line=inline('1/(4*x^2*(1 - 3/(4*x))^(2/3))');
[re1,step,x0_list]=NewtonDiedai(fx_line,dfx_line,0.78);
function [result,step,x0_list] = NewtonDiedai(fx,dfx,x0)
%此方法为牛顿迭代法求解非线性方程的根
% fx为原函数 dfx为导函数 x0为初值
% result为结果，step为迭代步数 x0_list为迭代结果的列表
x1=x0;
step=1;
x0_list=0;
while 1
x0=x1;
x1=x0-feval(fx,x0)/feval(dfx,x0); %feval()函数执行指定的函数
if abs(x0-x1)<0.0001 %判断迭代结果精度是否达到要求
break;
end
if(step>55) %当步数大于50时认为迭代不收敛，退出
break;
end
step=step+1; %记录步数
x0_list=[x0_list x0];
end
result=x0;
end
题目三
%设置矩阵维数分别为100 1000 10000
n=[100,1000,10000];
stepj_list=0; %存储Jacobi方法的迭代步数
stepgs_list=0; %存储Gauss-Seidel方法的迭代步数
for k=1:3
%制造A矩阵
a=zeros(n(k));
a(1,2)=-1;
a(n(k),n(k)-1)=-1;
for i=1:n(k)
a(i,i)=3;
end
for i=2:n(k)-1
a(i,i-1)=-1;
a(i,i+1)=-1;
end
%制造b向量
b=ones(n(k),1);
b(1)=2;
b(n(k))=2;
%分别调用函数，求解各自的方程解，以及迭代次数
[xj,countj]=Jacobi(a,b,zeros(n(k),1),10^(-6));
[xgs,countgs]=Gauss_Seidel(a,b,zeros(n(k),1),10^(-6));
stepj_list=[stepj_list,countj];
stepgs_list=[stepgs_list,countgs];
end
function [x,count] = Jacobi(A,b,x0,eps)
%利用jacobi计算Ax=b的解
% x为计算所得解
% count为迭代次数
% A，b分别为Ax=b的系数与常数 x0为初值，eps为精度
%分离D L U
D=diag(diag(A)); %A的对角矩阵 只有一个diag是一个列向量
L=-tril(A,-1); %利用tril求解下三角矩阵，-1代表不包含对角线
U=-triu(A,1); %利用triu求解上三角矩阵，1代表不包含对角线
M=D\(L+U); %计算迭代矩阵
B=D\b;
x=M*x0+B;
count=1;
while norm(x-x0)>=eps %循环，直至数值解满足精度要求
x0=x;
x=M*x0+B;
count=count+1;
if count>=10000 %当迭代次数过多，则不收敛
disp('Warning:迭代次数太多,仍未达精度要求,可能不收敛!')
return;
end
end
function [x,count] = Gauss_Seidel(A,b,x0,eps)
%GAUSS_SEIDEL 计算Ax=b的解
% x为计算所得解
% count为迭代次数
% A，b分别为Ax=b的系数与常数 x0为初值，eps为精度
%分离D L U
D=diag(diag(A)); %A的对角矩阵 只有一个diag是一个列向量
L=-tril(A,-1); %利用tril求解下三角矩阵，-1代表不包含对角线
U=-triu(A,1); %利用triu求解上三角矩阵，1代表不包含对角线
M=-(D+L)\U; %计算迭代矩阵
B=(D+L)\b;
x=M*x0+B;
count=1;
while norm(x-x0)>=eps %循环，直至数值解满足精度要求
x0=x;
x=M*x0+B;
count=count+1;
if count>=10000 %当迭代次数过多，则不收敛
disp('Warning:迭代次数太多,仍未达精度要求,可能不收敛!')
return;
end
end
end
题目四
%% 提供原函数图像的点
x_yuan=linspace(-1,1,100);
y_yuan=1./(1.+12.*x_yuan.*x_yuan);
%% 分别用4次、8次、16次，20次插值多项式
dist=2; %区间长度
list_poly_count=[4,8,16,20]; %多项式的次数
% 任意一个次数构造
%四次
step=dist/list_poly_count(1);
x0=-1:step:1; %计算需多少个节点
y0=1./(1.+12.*x0.*x0);
[C,L]= Lagrange(x0,y0); %调用函数计算
resultx4=-1:0.0001:1;
resulty4=polyval(C,resultx4); %计算resultx4对应的插值多项式的值，用于画图
%8次
step=dist/list_poly_count(2);
x0=-1:step:1; %计算需多少个节点
y0=1./(1.+12.*x0.*x0);
[C,L]= Lagrange(x0,y0); %调用函数计算
resultx8=-1:0.0001:1;
resulty8=polyval(C,resultx8); %计算resultx8对应的插值多项式的值，用于画图
%16次
step=dist/list_poly_count(3);
x0=-1:step:1; %计算需多少个节点
y0=1./(1.+12.*x0.*x0);
[C,L]= Lagrange(x0,y0);%调用函数计算
resultx16=-1:0.0001:1;
resulty16=polyval(C,resultx16);%计算resultx对应的插值多项式的值，用于画图
%20次
step=dist/list_poly_count(4);
x0=-1:step:1;
y0=1./(1.+12.*x0.*x0);%计算需多少个节点
[C,L]= Lagrange(x0,y0);%调用函数计算
resultx20=-1:0.0001:1;
resulty20=polyval(C,resultx20);%计算resultx对应的插值多项式的值，用于画图
%% 图像绘制
subplot(2,2,1);
plot(resultx4,resulty4,x_yuan,y_yuan,'linewidth',5);
title("\fontsize{20}4次插值与原函数对比图");
legend("4次插值","原函数图像");
subplot(2,2,2);
plot(resultx8,resulty8,x_yuan,y_yuan,'linewidth',5);
title("\fontsize{20}8次插值与原函数对比图");
legend("8次插值","原函数图像");
subplot(2,2,3);
plot(resultx16,resulty16,x_yuan,y_yuan,'linewidth',5);
title("\fontsize{20}16次插值与原函数对比图");
legend("16次插值","原函数图像");
subplot(2,2,4);
plot(resultx20,resulty20,x_yuan,y_yuan,'linewidth',5);
title("\fontsize{20}20次插值与原函数对比图");
legend("20次插值","原函数图像");
function [C,L]= Lagrange(X,Y)
% Lagrange用于求解拉格朗日插值
% X为各插值点的x坐标
% Y为各插值点的y坐标
% C为返回的插值多项式系数
% L为插值多项式（sym）
m = length(X); %求解共有多少插值点
for k = 1 : m
V = 1;
for i = 1 : m
if k ~= i
%求解每个基函数的每一项
V = conv(V,poly(X(i))) / (X(k) - X(i)); %用conv进行多项式的计算
%其中poly函数返回一个类似于x-xi的形式
end
end
L1(k, :) = V; %存储每个基函数的各个系数
l(k, :) = poly2sym(V);%将基函数的系数转为多项式形式，并存储
end
%分别计算多项式系数以及sym形式的多项式，方便展示
C = Y * L1;
L = Y * l;
end
根据图4-1，我们可以发现对于拉格朗日插值多项式，插值的次数越高，中间部分对原函数的贴合程度越高，但是Runge现象会非常明显；次数较低时中间部分的逼近效果会比较差，但是Runge现象也会减弱。
题目五
%% 传入数据
year={'1994','1995','1996','1997','1998','1999','2000','2001','2002','2003'};
year0=1994;
year_add=[0,1,2,3,4,5,6,7,8,9];
bbl=[67.052,68.008,69.803,72.024,73.400,72.063,74.669,74.487,74.065,76.777];
%% 拉格朗日插值
[C,L]= Lagrange(year_add,bbl);%为了方便计算，特将x化到0-9
predit_x=0:0.01:9;%计算需多少个节点
result=polyval(C,predit_x);%计算predit_x对应的插值多项式的值，用于画图
%绘制插值多项式图
hold on
plot(predit_x,result,'linewidth',5);
plot(year_add,bbl,"*",'linewidth',5);
set(gca,'XTicklabel',year);
title('\fontsize{20}9次插值多项式');
%% 求解2010年
add=2010-year0;
result_2010=polyval(C,add);%计算2010年对应的插值多项式的值
% 可以看出有Runge现象
%% 求解三次插值
year_1_4=year_add(1:4);
bbl_1_4=bbl(1:4);
%利用上述方法继续求解三次插值
[C_4,L_4,L1_4,l_4]= Lagrange(year_1_4,bbl_1_4);
predit_x=0:0.01:9;
result_1_4=polyval(C_4,predit_x);
hold on
plot(predit_x,result_1_4,'linewidth',5);
plot(year_add,bbl,"*",'linewidth',5);
set(gca,'XTicklabel',year);
title('\fontsize{20}3次插值多项式');
%% 用3次求1998年
add_98=1998-year0;
result_98=polyval(C_4,add_98);
% 没有Runge现象
%% 数据传入
year0=1994;
year_add=[0,1,2,3,4,5,6,7,8,9];
bbl=[67.052,68.008,69.803,72.024,73.400,72.063,74.669,74.487,74.065,76.777];
year={'1994','1995','1996','1997','1998','1999','2000','2001','2002','2003'};
%% 最小二乘直线拟合
a(:,1)=year_add';
a(:,2)=1;
%运用最小二乘原理计算
a_predit_=(a'*a)\(a'*bbl');
x_predit_=0:0.0001:9;
y_predit_=x_predit_*a_predit_(1)+a_predit_(2);
subplot(2,2,1);
hold on;
plot(x_predit_,y_predit_,'linewidth',5);
plot(year_add,bbl,"*",'linewidth',5);
set(gca,'XTicklabel',year);
title('\fontsize{20}最小二乘直线拟合');
%% 抛物线拟合
%运用polyfit进行二次拟合
c_predit_2=polyfit(year_add,bbl,2);
x_predit_2=0:0.0001:9;
%此处运用polyval计算拟合函数中各点坐标，用于绘图
y_predit_2=polyval(c_predit_2,x_predit_2);
subplot(2,2,2);
hold on;
plot(x_predit_2,y_predit_2,'linewidth',5);
plot(year_add,bbl,"*",'linewidth',5);
set(gca,'XTicklabel',year);
title('\fontsize{20}抛物线拟合');
%% 三次多项式拟合
%运用polyfit进行三次拟合
c_predit_3=polyfit(year_add,bbl,3);
x_predit_3=0:0.0001:9;
%此处运用polyval计算拟合函数中各点坐标，用于绘图
y_predit_3=polyval(c_predit_3,x_predit_3);
subplot(2,2,3);
hold on;
plot(x_predit_3,y_predit_3,'linewidth',5);
plot(year_add,bbl,"*",'linewidth',5);
set(gca,'XTicklabel',year);
title('\fontsize{20}三次多项式拟合');
%% 综合图
subplot(2,2,4);
hold on;
plot(x_predit_,y_predit_,x_predit_2,y_predit_2,x_predit_3,y_predit_3,'linewidth',2);
plot(year_add,bbl,"*",'linewidth',2);
legend('最小二乘','抛物线','三次多项式','原图像');
title('\fontsize{20}综合图像');
%% 预计2010年原油产量
add=2010-1994;
% 用最小二乘
predit_10_1=add*a_predit_(1)+a_predit_(2);
% 用抛物线拟合
predit_10_2=polyval(c_predit_2,add);
% 三次多项式拟合
predit_10_3=polyval(c_predit_3,add);
function [C,L]= Lagrange(X,Y)
% Lagrange用于求解拉格朗日插值
% X为各插值点的x坐标
% Y为各插值点的y坐标
% C为返回的插值多项式系数
% L为插值多项式（sym）
m = length(X); %求解共有多少插值点
for k = 1 : m
V = 1;
for i = 1 : m
if k ~= i
%求解每个基函数的每一项
V = conv(V,poly(X(i))) / (X(k) - X(i)); %用conv进行多项式的计算
%其中poly函数返回一个类似于x-xi的形式
end
end
L1(k, :) = V; %存储每个基函数的各个系数
l(k, :) = poly2sym(V);%将基函数的系数转为多项式形式，并存储
end
%分别计算多项式系数以及sym形式的多项式，方便展示
C = Y * L1;
L = Y * l;
end
经过分析图5-3以及表5-2，我们首先可以看出相比于插值多项式，三种类型的拟合多项式都可以比较好地预测原油产量。
通过分析图5-4，我们可以发现,虽然最小二乘直线能比较好地刻画原油产量随时间的变化，但是通过散点图我们可以看出点的分布在之后有弯的趋势，这更加符合二次或者三次多项式的特征；
在比较抛物线拟合与三次多项式拟合的图像后，我们可以发现三次多项式能顾及更多结点，而且走势也基本符合散点的走势，而且也没有出现过拟合的情况。
综上所示，我们可以得出三次多项式拟合曲线更能最好的表示表格中的数据。
题目六
%% 用复合梯形公式计算第一题
ft1=@(x)x^2*log(x);
[a,nn]=fuhuatixing(ft1,1,3,0.5*10^(-8));
%% 用复合梯形公式计算第二题
ft2=@(x)log(cos(x)+sin(x));
[a1,nn1]=fuhuatixing(ft2,0,pi/2,0.5*10^(-8));
%% 用复合梯形公式计算第三题
ft3=@(x)x/(2*exp(x)-exp(-x));
[a2,nn2]=fuhuatixing(ft3,0,1,0.5*10^(-8));
%% 用复合Simpson公式计算第一题
syms x;
f1=x^2*log(x);
[s,n] = Simpson(1,3,f1);
%% 用复合Simpson公式计算第二题
syms x;
f2=log(cos(x)+sin(x));
[s2,n2] = Simpson(0,pi/2,f2);
%% 用复合Simpson公式计算第三题
syms x;
f3=x/(2*exp(x)-exp(-x));
[s3,n3] = Simpson(0,1,f3);
%% 用Romberg公式计算第一题
[i1,num1]=Romberg(1,3,f1);
%% 用Romberg公式计算第二题
[i2,num2]=Romberg(0,pi/2,f2);
%% 用Romberg公式计算第三题
[i3,num3]=Romberg(0,1,f3);
function [A,N] = fuhuatixing(f1,a,b,q)
%用复化梯形公式进行积分求解
% a b：积分上下限 q:精度
%f1: 为所求函数的句柄
%A：积分结果
%N：小区间数量
m=1;
h=b-a;
T(1)=h/2*(f1(a)+f1(b));%初始值设定
%循环划分区间
while 1
h=h/2;
sum=0;
for j=1:2^(m-1)
sum=sum+f1(a+(2*j-1)*h);
end
S(1)=T(1)/2+h*sum;
for j=1:m
S(j+1)=S(j)+(S(j)-T(j))/(4^j-1);
end
if(abs(S(m+1)-T(m))<=q) %直至精度满足要求
A=S(m+1);
break;
end
T=S;
m=m+1;
end
%计算区间数量
N=(b-a)/h;
end
function [i,num] = Romberg(low,up,L)
%利用ROMBERG 方法进行复化求积
% low,up是积分上下限
% L是公式（sym的形式）
% i为积分结果，num为区间数量
tolerance=0.5*10^(-8); %限差
h=(up-low)/2;
t(1,1)=h*(eval(subs(L,low))+eval(subs(L,up)));
%subs()函数 表示将符号表达式中的某些符号变量替换为指定的新的变量
k=1;n=1;
m=1;
%根据书中的伪代码，可以写出复化函数
while 1
f=0;
for i=1:n
f=f+eval(subs(L,(low+(2*i-1)*h)));
end
t(k+1,1)=t(k,1)/2+h*f;
for m=1:k
t(k-m+1,m+1)=(4^m*t(k-m+2,m)-t(k-m+1,m))/(4^m-1);
end
if abs(t(1,m+1)-t(1,m))<tolerance
break;
end
%将区间分半，继续求解
h=h/2;
n=2*n;
k=k+1;
end
%取出romberg计算矩阵的最外侧数值即为最终积分解
i=t(1,m);
%计算区间数量
num=(up-low)/h;
function [ss,num] = Simpson(low,up,L)
% Simpson方法进行复化求积
% low,up是积分上下限
% L是公式（sym的形式）
% ss为积分结果，num为区间数量
tolerance=0.5*10^(-8); %限差
%进行初值计算
f1=eval(subs(L,low)+subs(L,up));
f2=eval(subs(L,(low+up)/2));
s0=(up-low)*(f1+4*f2)/6;
n=2;
h=(up-low)/4;
st=0;
%根据书中的伪代码，可以写出复化函数
while(1)
f3=0;
for i=1:n
f3=f3+eval(subs(L,low+(2*i-1)*h));
end
s=h*(f1+2*f2+4*f3)/3;
if abs(s-s0)<15*tolerance %直至精度满足要求
ss=s;
break;
end
st=st+1;
%将区间分半，继续求解
n=2*n;
h=h/2;
f2=f2+f3;
s0=s;
end
%计算区间数量
num=(up-low)/h;
end
根据表6-1、6-2、6-3所示，我们可以看出，两个方法在求不同的积分时，最终求得的值相等。
但Romberg求积方法是在复化梯形求积公式的基础上，应用了Richardson外推法构造的一种算法，是一种加速计算方法，故而与复合梯形公式/Simpson公式相比，所需的区间数量更少，能更快地求解。
题目七
%首先用dsolve('Dy=x','x')求解准确函数
%% 准确函数绘制
x0=0:0.001:1;
y0=x0.^2./2+1;
%% 用euler法求解第一个函数
%根据公式依次求解被积函数的积分值
% 分别按照步长0.1 0.05 0.025计算求解
% 步长为0.1求解
step=0.1;
x1=0:step:1;
y1=zeros(1,1/step+1);%存储每次迭代值
y1(1)=1;%初值为1
for i=2:size(y1,2) %euler法只与前一项有关，循环求解即可
y1(i)=y1(i-1)+step*x1(i-1);
end
% 0.05
step=0.05;
x2=0:step:1;
y2=zeros(1,1/step+1);
y2(1)=1;
for i=2:size(y2,2)
y2(i)=y2(i-1)+step*x2(i-1);
end
% 0.025
step=0.025;
x3=0:step:1;
y3=zeros(1,1/step+1);
y3(1)=1;
for i=2:size(y3,2)
y3(i)=y3(i-1)+step*x3(i-1);
end
%% 绘制用Euler法求解y'=x图
hold on;
plot(x0,y0,'linewidth',2);
plot(x1,y1,'linewidth',2);
plot(x2,y2,'linewidth',2);
plot(x3,y3,'linewidth',2);
legend('原函数','步长为0.1','步长为0.05','步长为0.025','Location','Best' );
title("\fontsize{20}用Euler法求解y'=x");
%% 用改进的euler法求解
%根据公式依次求解被积函数的积分值
% 分别按照步长0.1 0.05 0.025计算求解
% 步长为0.1求解
step=0.1;
x1=0:step:1;
y1=zeros(1,1/step+1);%存储每次迭代值
y1(1)=1;%初值为1
y11=y1;%存储y''
for i=2:size(y1,2)
y1(i)=y1(i-1)+step*x1(i-1);
y11(i)=y11(i-1)+step/2*(x1(i-1)+x1(i));
end
% 0.05
step=0.05;
x2=0:step:1;
y2=zeros(1,1/step+1);
y2(1)=1;
y21=y2;
for i=2:size(y2,2)
y2(i)=y2(i-1)+step*x2(i-1);
%首先运用euler做预估校正
y21(i)=y21(i-1)+step/2*(x2(i-1)+x2(i));
%将预估校正的值代入进行迭代，所得为改进的euler法所求的积分值
end
% 0.025
step=0.025;
x3=0:step:1;
y3=zeros(1,1/step+1);
y3(1)=1;
y31=y3;
for i=2:size(y3,2)
y3(i)=y3(i-1)+step*x3(i-1);
y31(i)=y31(i-1)+step/2*(x3(i-1)+x3(i));
end
%% 绘制用改进的Euler法求解y'=x图
hold on;
plot(x0,y0,'linewidth',2);
plot(x1,y1,'linewidth',2);
plot(x2,y2,'linewidth',2);
plot(x3,y3,'linewidth',2);
legend('原函数','步长为0.1','步长为0.05','步长为0.025','Location','Best' );
title("\fontsize{20}用改进的Euler法求解y'=x");
%首先用dsolve('Dy=x^2*y','x')求解准确函数
%% 准确函数
x0=0:0.001:1;
y0=exp(x0.^3./3);
%% 用euler法求解
%根据公式依次求解被积函数的积分值
% 分别按照步长0.1 0.05 0.025计算求解
% 0.1
step=0.1;
x1=0:step:1;
y1=zeros(1,1/step+1);%存储每次迭代值
y1(1)=1;%初值为1
for i=2:size(y1,2)
y1(i)=y1(i-1)+step*x1(i-1)^2*y1(i-1);
end%euler法只与前一项有关，循环求解即可
% 0.05
step=0.05;
x2=0:step:1;
y2=zeros(1,1/step+1);
y2(1)=1;
for i=2:size(y2,2)
y2(i)=y2(i-1)+step*x2(i-1)^2*y2(i-1);
end
% 0.025
step=0.025;
x3=0:step:1;
y3=zeros(1,1/step+1);
y3(1)=1;
for i=2:size(y3,2)
y3(i)=y3(i-1)+step*x3(i-1)^2*y3(i-1);
end
%% 绘制用Euler法求解y'=x^{2}y图
hold on;
plot(x0,y0,'linewidth',2);
plot(x1,y1,'linewidth',2);
plot(x2,y2,'linewidth',2);
plot(x3,y3,'linewidth',2);
legend('原函数','步长为0.1','步长为0.05','步长为0.025','Location','Best' );
title("\fontsize{20}用Euler法求解y'=x^{2}y");
%% 用改进的euler法求解
% 步长0.1 0.05 0.025
% 0.1
step=0.1;
x1=0:step:1;
y1=zeros(1,1/step+1);
y1(1)=1;
y11=y1;
for i=2:size(y1,2)
y1(i)=y1(i-1)+step*x1(i-1)^2*y1(i-1);%首先运用euler做预估校正
y11(i)=y11(i-1)+step/2*(x1(i-1)^2*y11(i-1)+x1(i)^2*y1(i));
%将预估校正的值代入进行迭代，所得为改进的euler法所求的积分值
end
% 0.05
step=0.05;
x2=0:step:1;
y2=zeros(1,1/step+1);%存储每次迭代值
y2(1)=1;%初值为1
y21=y2;%存储y''
for i=2:size(y2,2)
y2(i)=y2(i-1)+step*x2(i-1)^2*y2(i-1);
y21(i)=y21(i-1)+step/2*(x2(i-1)^2*y21(i-1)+x2(i)^2*y2(i));
end
% 0.025
step=0.025;
x3=0:step:1;
y3=zeros(1,1/step+1);
y3(1)=1;
y31=y3;
for i=2:size(y3,2)
y3(i)=y3(i-1)+step*x3(i-1)^2*y3(i-1);
y31(i)=y31(i-1)+step/2*(x3(i-1)^2*y31(i-1)+x3(i)^2*y3(i));
end
%% 绘制用改进的Euler法求解y'=x^{2}y图
hold on;
plot(x0,y0,'linewidth',2);
plot(x1,y1,'linewidth',2);
plot(x2,y2,'linewidth',2);
plot(x3,y3,'linewidth',2);
legend('原函数','步长为0.1','步长为0.05','步长为0.025','Location','Best' );
title("\fontsize{20}用改进的Euler法求解y'=x^{2}y");
%首先用dsolve('Dy=2*(x+1)*y','x')求解准确函数
%% 准确函数
x0=0:0.001:1;
y0=exp(x0.*(x0+ 2));
%% 用euler法求解
%根据公式依次求解被积函数的积分值
% 分别按照步长0.1 0.05 0.025计算求解
%步长： 0.1
step=0.1;
x1=0:step:1;
y1=zeros(1,1/step+1);%存储每次迭代值
y1(1)=1;%初值为1
for i=2:size(y1,2)
y1(i)=y1(i-1)+step*2*(x1(i-1)+1)*y1(i-1);
end
% 0.05
step=0.05;
x2=0:step:1;
y2=zeros(1,1/step+1);
y2(1)=1;
for i=2:size(y2,2)%euler法只与前一项有关，循环求解即可
y2(i)=y2(i-1)+step*2*(x2(i-1)+1)*y2(i-1);
end
% 0.025
step=0.025;
x3=0:step:1;
y3=zeros(1,1/step+1);
y3(1)=1;
for i=2:size(y3,2)
y3(i)=y3(i-1)+step*2*(x3(i-1)+1)*y3(i-1);
end
%% 绘制用Euler法求解y'=2(x+1)y图
hold on;
plot(x0,y0,'linewidth',2);
plot(x1,y1,'linewidth',2);
plot(x2,y2,'linewidth',2);
plot(x3,y3,'linewidth',2);
legend('原函数','步长为0.1','步长为0.05','步长为0.025','Location','Best' );
title("\fontsize{20}用Euler法求解y'=2(x+1)y");
%% 用改进的euler法求解
% 步长0.1 0.05 0.025
% 0.1
step=0.1;
x1=0:step:1;
y1=zeros(1,1/step+1);
y1(1)=1;
y11=y1;
for i=2:size(y1,2)
y1(i)=y1(i-1)+step*2*(x1(i-1)+1)*y1(i-1);%首先运用euler做预估校正
y11(i)=y11(i-1)+step/2*(2*y11(i-1)*(x1(i-1)+1)+2*y1(i)*(x1(i)+1));
%将预估校正的值代入进行迭代，所得为改进的euler法所求的积分值
end
% 0.05
step=0.05;
x2=0:step:1;
y2=zeros(1,1/step+1);
y2(1)=1;
y21=y2;
for i=2:size(y2,2)
y2(i)=y2(i-1)+step*2*(x2(i-1)+1)*y2(i-1);
y21(i)=y2(i-1)+step/2*(2*y21(i-1)*(x2(i-1)+1)+2*y2(i)*(x2(i)+1));
end
% 0.025
step=0.025;
x3=0:step:1;
y3=zeros(1,1/step+1);%存储每次迭代值
y3(1)=1;%初值为1
y31=y3;%存储y''
for i=2:size(y3,2)
y3(i)=y3(i-1)+step*2*(x3(i-1)+1)*y3(i-1);
y31(i)=y3(i-1)+step/2*(2*y31(i-1)*(x3(i-1)+1)+2*y3(i)*(x3(i)+1));
end
%% 绘制用改进的Euler法求解y'=2(x+1)y图
hold on;
plot(x0,y0,'linewidth',2);
plot(x1,y1,'linewidth',2);
plot(x2,y2,'linewidth',2);
plot(x3,y3,'linewidth',2);
legend('原函数','步长为0.1','步长为0.05','步长为0.025','Location','Best' );
title("\fontsize{20}用改进的Euler法求解y'=2(x+1)y");
6131
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
