一.数值积分
常用解法：Newton-Cotes型求积公式。当节点为两个时为梯形公式，节点为三个时为Simpson公式
使用高阶NC公式时虽然精度较高，但是稳定性得不到保证，低阶则相反。为此我们将积分区间划分成等长的子区间，在每个子区间上使用低阶的NC公式进行计算，最后再求和得到近似值
复化梯形公式：我们知道，将后式的每一个积分值以中间值乘以区间长度代替：
便得到复化梯形公式
复化Simpson公式：在每个子区间上利用Simpson公式就有
下面是一个简单的例子：我们通过公式对
的值进行估算
syms x
f=@(x)[4/(1+x^2)];
for i=1:10
n=2^i;
h(1,i)=1/n;%h为步长，这里我们从0到1积分
x=linspace(0,1,n+1);%注意这里的n+1是节点数
Tn=0;
Sn=0;
error1(1,i)=0;
error2(1,i)=0;
for k=1:n
Tn=Tn+(h(1,i)/2)*(f(x(1,k))+f(x(1,k+1)));
Sn=Sn+(h(1,i)/6)*(f(x(1,k))+4*f((x(1,k)+x(1,k+1))/2)+f(x(1,k+1)));
end
fprintf('第%d次复化梯形公式计算结果为：%.10f\n',i,Tn);
error1(1,i)=abs(Tn-pi);
fprintf('误差为：%.10f\n',error1(1,i));
fprintf('第%d次复化Simpson公式计算结果为：%.10f\n',i,Sn);
error2(1,i)=abs(Sn-pi);
fprintf('误差为：%.10f\n',error2(1,i));
end
Romberg(0,1,f);
figure(1),plot(log(h),log(error1));
figure(2),plot(log(h),log(error2));
区间长度每进行一次便缩短一半，通过不断改变h的值可以得到两种方法误差与h之间的对数关系图
理论上来说收敛阶分别为2阶和4阶，但第二个收敛速度明显快于4阶，可能是由于函数选取较特殊。同时，继续缩短步长h后出现了误差回弹的情况，这是因为我们在子区间上使用了低阶的Newton-Cotes公式来拟合积分值，当区间越来越多达到一个界限时，截断误差相比理论上精度的提高会起到主要影响（理论上来说越小越好），并且matlab中使用double型小数，并不能保证高精度的准确性。
Romberg算法：在进行过梯形的二等分求积之后再对每个区间各二等分一次，步长h不断缩短（在原来的分割基础上不断细化）递推关系式为
龙贝格算法又被称作变步长方法，我们直接使用矩阵来递推求解即可。需要注意的是，为了加速收敛，我们使用了Richardson外推加速法。精度设置为1e-10
function Romberg(a,b,f)
T(1,1)=(b-a)/2*(f(a)+f(b));
k=1;
epsilon=1e-10;
error3=1;
while(error3>epsilon)
sum=0;
for i=1:2^(k-1)
sum=sum+f(a+(2*i-1)*(b-a)/2^k);
end
T(1,k+1)=0.5*T(1,k)+((b-a)/2^k)*sum;
for m=1:k
T(m+1,k+1)=4^m/(4^m-1)*T(m,k+1)-1/(4^m-1)*T(m,k);
end
k=k+1;
error3=abs(T(k,k)-T(k-1,k-1));
end
fprintf('Romberg进行了%d次二等分,近似结果为%.10f，误差为%.11f',k-1,T(k,k),error3);
二.常微分方程初值问题
对于常规一阶常微分方程组
我们考虑以Taylor展开为核心的Euler方法，它将一阶导数通过方程转化为
得到近似计算公式
通过隐式方程改进的Euler公式
将斜率值加权组合，我们得到了更高精度的Runge-Kutta方法。根据参数个数的选取，可以得到不同阶的R-K迭代方程组，综合计算量与精度方面的考虑，我们一般使用四阶R-K公式：
来看一个具体例子
显式欧拉法
function Euler1(a,b,h)
n=(b-a)/h;
y(1,1)=0;
x(1,1)=a;
for k=1:n-1
x(1,k+1)=a+k*h;
y(1,k+1)=y(1,k)+h*g(x(1,k),y(1,k));
end
disp('近似解:');
disp(y);
ry=x.^(-2).*(exp(x)-exp(1));
disp('误差：');
disp(abs(y-ry));
plot(x,y,'r',x,ry,'g');
改进欧拉法
function Euler2(a,b,h)
n=(b-a)/h;
y(1,1)=0;
x(1,1)=a;
for k=1:n-1
x(1,k+1)=a+k*h;
yba(1,k+1)=y(1,k)+h*g(x(1,k),y(1,k));
y(1,k+1)=y(1,k)+h/2*(g(x(1,k),y(1,k))+g(x(1,k+1),yba(1,k+1)));
end
disp('近似解:');
disp(y);
ry=x.^(-2).*(exp(x)-exp(1));
disp('误差：');
disp(abs(y-ry));
plot(x,y,'r',x,ry,'g')；
R-K方法
function RK4(a,b,h)
n=(b-a)/h;
y(1,1)=0;
x(1,1)=a;
for k=1:n-1
x(1,k+1)=a+k*h;
k1=g(x(1,k),y(1,k));
k2=g(x(1,k)+h/2,y(1,k)+h/2*k1);
k3=g(x(1,k)+h/2,y(1,k)+h/2*k2);
k4=g(x(1,k)+h,y(1,k)+h*k3);
y(1,k+1)=y(1,k)+h/6*(k1+2*k2+2*k3+k4);
end
disp('近似解:');
disp(y);
ry=x.^(-2).*(exp(x)-exp(1));
disp('误差：');
disp(abs(y-ry));
plot(x,y,'r',x,ry,'g');
其中前两个方法均为1阶精度，RK为4阶精度，RK在这个例子中的精度也是最高的
经过对1<x<b中b的取值调整（固定步长为0.1），b在<=20时精度都压缩在1e-2内，但b=30时靠近b的x拟合的误差便达到了三位数，b=35时误差为1e+4级，而b=40后误差达到了1e+6级
分析图像我们猜测这可能是y的取值过大造成的误差，因为截断误差允许的精度范围并不大
查阅资料后发现理论来讲RK算法存在稳定域
由此我们知道确实存在这样一个临界步长使得算法不稳定。
