Running Snail

function [x, y] = improvedEuler(fname, inte, y0, h)x=inte(1):h:inte(2);y(1)=y0;for n=1:length(x)-1 y(n)=vpa(y(n),6); k1=feval(fname,x(n),y(n)); y(n+1)=y(n)+h*k1; k2=feval(fname,x(n+1),y(n+1)); y(n+1)=y(n)+h*(k1+k2)/2;endx=x'y=y'..

用变步长的中点方法和三点求导公式求exe^xex在x=1初的导数值，并比较步长的变化对解的影响function T = var_size(fname,h)format longk=10;for i=0:k m=h/(2^i); T=feval(fname,m); disp([i,T])endfname=inline('(exp(1+x)-exp(1-x))/(2*x)');h0=1;h1=0.1;h2=0.01;var_size(fname,h0)var_size(

利用复化辛普森求积公式计算∫abf(x)dx\int _ { a } ^ { b } f ( x ) d x∫ab​f(x)dx的近似值辛普森求积function s=simpson( f_name,a, b, n)%复化辛普生公式% f_name为要求的定函数y=f(x)所在的程序文件名% a为积分下限% b为积分上限% n为积分区间[a,b]划分成小区间的等份数h=(b-a)/n;fa=feval(f_name,a);fb=feval(f_name,b);s=fa-fb;x=a

function x = agui_gauss(a,b)n=length(b);a=[a,b];for k=1:(n-1) [ar,r]=max(abs(a(k:n,k))); r=r+k-1; if r>k t=a(k,:);a(k,:)=a(r,:);a(r,:)=t; end % 消元 a((k+1):n,(k+1):(n+1)) = a((k+1):n,(k+1):(n+1))-a((k+1):n,k)/a(k,k)*a(.

agui_lagrange.m：function f=agui_lagrange(x0,y0,x)% x0为节点向量，y0为节点上的函数值，x为插值点，f为返回插值n=length(x0);m=length(x);format longs=0.0;for k=1:n p=1.0; for j=1:n if j~=k p=p*(x-x0(j))/(x0(k)-x0(j)); end end s=p*y0(k).

Agui_GS.m：function x=Agui_GS(A,b)%方程Ax=b，x0为初始向量%e为精度，N为最大迭代次数n=length(b);N=100;e=1e-4;x0=zeros(n,1);x=x0;x0=x+2*e;k=0;A1=tril(A);A2=inv(A1);while norm(x0-x,inf)>e&&k<N k=k+1; x0=x; x=-A2*(A-A1)*x0+A2*b; format l.

function [x,y]=after_method()%数组a存储三角矩阵A的主对角线元素，c、d存储主对角线上边下边带宽为1的元素%追赶法 a=[2,2,2,2,2]; c=[-1,-1,-1,-1]; d=[-1,-1,-1,-1]; b=[1,0,0,0,0]; n=length(a); n1=length(c); n2=length(d); %错误检查 if n1~=n2%存储矩阵的数组维数错误 err.

function [l,u,x,y] = agui_lu(a,b)% 求可逆矩阵的分解，l为下三角矩阵，u为上三角矩阵n = length(a);u = zeros(n,n);l = eye(n,n);u(1,:)=a(1,:);l(2:n,1)=a(2:n,1)/u(1,1);for k=2:n u(k,k:n)=a(k,k:n)-l(k,1:k-1)*u(l:k-1,k:n); l(k+1:n,k)=( a(k+1:n,k)-l(k+1:n,1:k-1)*u(1:k-1,k

MATLAB函数文件function x=Agui_GS(A,b)%方程Ax=b，x0为初始向量%e为精度，N为最大迭代次数n=length(b);N=100;e=1e-4;x0=zeros(n,1);x=x0;x0=x+2*e;k=0;A1=tril(A);A2=inv(A1);while norm(x0-x,inf)>e&&k<N k=k+1; x0=x; x=-A2*(A-A1)*x0+A2*b; format lo

Matlab代码函数文件function x = agui_jacobi(a,b)%ax=b，a为系数矩阵，b为右端向量，x0为初始向量n=length(b); % 精度N=100; % 最大迭代次数e=1e-6;x0=zeros(n,1); % 初始向量默认为0向量x=x0;x0=x+2*e;k=0; D=diag(diag(a)); % 调用 diag 两次将返回一个包含a的对角元素的对角矩阵

可以看到刚开始小数只有四位修改为long就可以了

function x = agui_newton(fname,dfname,x0,e)N=10;x=x0;x0=x+2*e;k=0;while abs(x0-x)>e&&k<N k=k+1; x0=x; x=x0-feval(fname,x0)/feval(dfname,x0); disp(x)endif k==N warning('已达最大迭代次数');endend命令行>> fun=inline

function x = agui_bisect(fname,a,b,e)fa = feval(fname,a);fb = feval(fname,b);if fa*fb>0 error("两端函数值同号");endk=0;x=(a+b)/2;while(b-a)>(2*e) fx=feval(fname,x); if fa*fx<0 b=x; fb=fx; else a=x; f