有限基础分析

复制以下各段代码，打开matlab脚本文件（.m文件），粘贴后直接运行即可。例2-1代码%二维插值示例clearall;closeall;clc;%定义二维数据点x=[1.2,2.1,3.5,4.8,5.2,6.7,7.3,8.9,9.1,10.0];%x坐标y=[2.5,3.7,1.8,4.2,5.9,6.1,7.8,8.3,9.5,10.2];%y坐标%创建参数t表示点的顺序t=1:length(x);%创建细粒度参数用于插值ti=linspace(min(t),max(t),100);legend('原始数据点','线性插值');%对x和y分别进行插值interpolation_method='linear';%线性插值xi=interp1(t,x,ti,interpolation_method);yi=interp1(t,y,ti,interpolation_method);%可视化原始数据点和插值路径figure('Position',[100,100,800,600]);%原始数据点和插值路径plot(x,y,'ro','MarkerSize',8,'MarkerFaceColor','r');holdon;plot(xi,yi,'b-','LineWidth',2);xlabel('X坐标');ylabel('Y坐标');gridon;holdoff;例2-2代码%二维插值示例clearall;closeall;clc;%定义二维数据点x=[1.2,2.1,3.5,4.8,5.2,6.7,7.3,8.9,9.1,10.0];%x坐标y=[2.5,3.7,1.8,4.2,5.9,6.1,7.8,8.3,9.5,10.2];%y坐标%创建参数t表示点的顺序t=1:length(x);%创建细粒度参数用于插值ti=linspace(min(t),max(t),100);%对x和y分别进行插值interpolation_method='cubic';%三次函数插值xi=interp1(t,x,ti,interpolation_method);yi=interp1(t,y,ti,interpolation_method);%可视化原始数据点和插值路径figure('Position',[100,100,800,600]);%原始数据点和插值路径plot(x,y,'ro','MarkerSize',8,'MarkerFaceColor','r');holdon;plot(xi,yi,'b-','LineWidth',2);xlabel('X坐标');ylabel('Y坐标');gridon;legend('原始数据点','三次函数插值');holdoff;例2-3代码%二维插值示例clearall;closeall;clc;%定义二维数据点x=[1.2,2.1,3.5,4.8,5.2,6.7,7.3,8.9,9.1,10.0];%x坐标y=[2.5,3.7,1.8,4.2,5.9,6.1,7.8,8.3,9.5,10.2];%y坐标%创建参数t表示点的顺序t=1:length(x);%创建细粒度参数用于插值gridon;ti=linspace(min(t),max(t),100);%对x和y分别进行插值interpolation_method='pchip';%三次函数分段插值xi=interp1(t,x,ti,interpolation_method);yi=interp1(t,y,ti,interpolation_method);%可视化原始数据点和插值路径figure('Position',[100,100,800,600]);%原始数据点和插值路径plot(x,y,'ro','MarkerSize',8,'MarkerFaceColor','r');holdon;plot(xi,yi,'b-','LineWidth',2);xlabel('X坐标');ylabel('Y坐标');legend('原始数据点','三次函数分段插值');holdoff;例2-4代码%二维数据线性拟合示例clearall;closeall;clc;%定义二维数据点（仅X和Y坐标）x=[1.2,2.1,3.5,4.8,5.2,6.7,7.3,8.9,9.1,10.0];%x坐标y=[2.5,3.7,1.8,4.2,5.9,6.1,7.8,8.3,9.5,10.2];%y坐标%拟合p=polyfit(x,y,1);%一次多项式（线性）拟合，数字越大，阶数越高y1=polyval(p,x);plot(x,y,'ro','MarkerSize',8,'MarkerFaceColor','r');holdon;plot(x,y1,'b-','LineWidth',2)holdoffxlabel('X坐标');ylabel('Y坐标');gridon;legend('原始数据点','线性拟合');例2-5代码%数值微分示例clearall;closeall;clc;%定义函数和采样点x=linspace(0,2*pi,100);%采样点y=sin(x);%原函数值%计算理论导数（解析解）dy_exact=cos(x);%sin(x)的导数是cos(x)%计算步长h=x(2)-x(1);%采样间隔%使用不同方法计算数值导数dy_forward=zeros(size(y));%向前差分dy_backward=zeros(size(y));%向后差分dy_central=zeros(size(y));%中心差分%向前差分(i=1到n-1)fori=1:length(y)-1dy_forward(i)=(y(i+1)-y(i))/h;enddy_forward(end)=dy_forward(end-1);%最后一点使用前一点的值%向后差分(i=2到n)fori=2:length(y)dy_backward(i)=(y(i)-y(i-1))/h;enddy_backward(1)=dy_backward(2);%第一点使用后一点的值%中心差分(i=2到n-1)fori=2:length(y)-1dy_central(i)=(y(i+1)-y(i-1))/(2*h);enddy_central(1)=dy_central(2);%第一点使用后一点的值dy_central(end)=dy_central(end-1);%最后一点使用前一点的值%计算误差error_forward=abs(dy_forward-dy_exact);error_backward=abs(dy_backward-dy_exact);error_central=abs(dy_central-dy_exact);%可视化结果figure('Position',[100,100,1200,800]);%绘制原函数subplot(2,2,1);plot(x,y,'b-','LineWidth',2);title('原函数y=sin(x)');xlabel('x');ylabel('y');gridon;%绘制导数subplot(2,2,2);plot(x,dy_exact,'k-','LineWidth',2,'DisplayName','解析解');holdon;plot(x,dy_forward,'r--','LineWidth',1.5,'DisplayName','向前差分');plot(x,dy_backward,'g--','LineWidth',1.5,'DisplayName','向后差分');plot(x,dy_central,'m--','LineWidth',1.5,'DisplayName','中心差分');title('导数计算结果对比');xlabel('x');ylabel('dy/dx');gridon;legend;%绘制误差subplot(2,2,3);semilogy(x,error_forward,'r-','LineWidth',1.5,'DisplayName','向前差分误差');holdon;semilogy(x,error_backward,'g-','LineWidth',1.5,'DisplayName','向后差分误差');semilogy(x,error_central,'m-','LineWidth',1.5,'DisplayName','中心差分误差');title('数值微分误差对比(对数刻度)');xlabel('x');ylabel('绝对误差');gridon;legend;%误差统计subplot(2,2,4);bar([123],[mean(error_forward)mean(error_backward)mean(error_central)]);title('平均误差比较');xticklabels({'向前差分','向后差分','中心差分'});ylabel('平均绝对误差');gridon;例2-6代码%离散数据点的积分示例clearall;closeall;clc;%定义被积函数（仅用于生成离散数据，实际应用中可能只有数据点）f=@(x)sin(x);%被积函数：f(x)=sin(x)a=0;%积分下限b=pi;%积分上限I_exact=2;%解析解%生成不同密度的离散数据点n_values=[5,10,20,50,100];%采样点数num_methods=3;%三种积分方法errors=zeros(length(n_values),num_methods);%对每个采样点数进行计算fori=1:length(n_values)n=n_values(i);x=linspace(a,b,n);%均匀分布的采样点y=f(x);%对应的函数值%1.矩形法（左端点）h=(b-a)/(n-1);%步长I_rect=h*sum(y(1:end-1));errors(i,1)=abs(I_rect-I_exact);%2.梯形法（MATLAB内置函数trapz）I_trap=trapz(x,y);errors(i,2)=abs(I_trap-I_exact);%3.辛普森法（要求n为奇数，即偶数个区间）ifmod(n,2)==1I_simp=simpson(x,y);%调用自定义Simpson函数else%如果n为偶数，去掉最后一个点使其变为奇数I_simp=simpson(x(1:end-1),y(1:end-1));enderrors(i,3)=abs(I_simp-I_exact);end%可视化结果figure('Position',[100,100,1200,800]);%绘制被积函数和离散点（以n=20为例）subplot(2,2,1);x_plot=linspace(a,b,1000);y_plot=f(x_plot);plot(x_plot,y_plot,'b-','LineWidth',2);holdon;n_example=20;x_example=linspace(a,b,n_example);y_example=f(x_example);scatter(x_example,y_example,50,'ro','filled');title(sprintf('被积函数与离散数据点(n=%d)',n_example));xlabel('x');ylabel('f(x)');gridon;legend('解析函数','离散数据点');%绘制梯形积分示例subplot(2,2,2);plot(x_plot,y_plot,'b-','LineWidth',2);holdon;forj=1:length(x_example)-1fill([x_example(j),x_example(j+1),x_example(j+1),x_example(j)],...[0,0,y_example(j+1),y_example(j)],'r','FaceAlpha',0.3);endtitle(sprintf('梯形积分近似(n=%d)',n_example));xlabel('x');ylabel('f(x)');gridon;legend('解析函数','梯形近似');%绘制误差比较subplot(2,2,3);loglog(n_values,errors(:,1),'ro-','LineWidth',1.5,'DisplayName','矩形法');holdonloglog(n_values,errors(:,2),'bo--','LineWidth',1.5,'DisplayName','梯形法');loglog(n_values,errors(:,3),'go-','LineWidth',1.5,'DisplayName','辛普森法');holdofftitle('离散数据积分误差比较');xlabel('采样点数n');ylabel('绝对误差');gridon;legend;%绘制误差收敛率subplot(2,2,4);loglog(n_values,errors(:,1),'ro-','LineWidth',1.5,'DisplayName','矩形法');loglog(n_values,errors(:,2),'bo-','LineWidth',1.5,'DisplayName','梯形法');loglog(n_values,errors(:,3),'go-','LineWidth',1.5,'DisplayName','辛普森法');%理论收敛率参考线ref_n=linspace(min(n_values),max(n_values),100);ref_h=(b-a)./(ref_n-1);loglog(ref_n,0.5*ref_h,'r--','LineWidth',1,'DisplayName','O(h)');holdonloglog(ref_n,0.1*ref_h.^2,'b--','LineWidth',1,'DisplayName','O(h²)');loglog(ref_n,0.01*ref_h.^4,'g--','LineWidth',1,'DisplayName','O(h⁴)');holdofftitle('误差收敛率对比');xlabel('采样点数n');ylabel('绝对误差');gridon;legend;axistight;%自定义Simpson函数（处理离散数据）functionI=simpson(x,y)n=length(x);ifmod(n,2)==0error('Simpson法要求数据点数n为奇数（偶数个区间）');endh=(x(end)-x(1))/(n-1);I=h/3*(y(1)+y(end)+4*sum(y(2:2:end-1))+2*sum(y(3:2:end-2)));end例2-7代码%线性方程组求解示例clearall;closeall;clc;%定义线性方程组Ax=bA=[4,-1,1;-1,4,-2;1,-2,4];%系数矩阵（对称正定）b=[12;-1;5];%常数向量%方法1：左除运算符（推荐方法）x1=A\b;%方法2：矩阵求逆（不推荐，计算效率低）x2=inv(A)*b;%方法3：LU分解[L,U]=lu(A);%LU分解y=L\b;%解Ly=bx3=U\y;%解Ux=y%方法4：Cholesky分解（适用于对称正定矩阵）R=chol(A);%Cholesky分解，A=R'*Ry=R'\b;%解R'y=bx4=R\y;%解Rx=y%方法5：QR分解[Q,R]=qr(A);%QR分解y=Q'*b;%解Qy=bx5=R\y;%解Rx=y%验证解的正确性err1=norm(A*x1-b);err2=norm(A*x2-b);err3=norm(A*x3-b);err4=norm(A*x4-b);err5=norm(A*x5-b);%显示结果fprintf('线性方程组求解结果：\n');fprintf('----------------------------------------\n');fprintf('方法1（左除）：x=[%.6f,%.6f,%.6f],误差=%.10f\n',x1,err1);fprintf('\n');fprintf('方法2（求逆）：x=[%.6f,%.6f,%.6f],误差=%.10f\n',x2,err2);fprintf('\n');fprintf('方法3（LU分解）：x=[%.6f,%.6f,%.6f],误差=%.10f\n',x3,err3);fprintf('\n');fprintf('方法4（Cholesky分解）：x=[%.6f,%.6f,%.6f],误差=%.10f\n',x4,err4);fprintf('\n');fprintf('方法5（QR分解）：x=[%.6f,%.6f,%.6f],误差=%.10f\n',x5,err5);fprintf('----------------------------------------\n');例2-8代码%线性方程组求解示例clearall;closeall;clc;%定义线性方程组Ax=bA=[4,-1,1;-1,4,-2;1,-2,4];%系数矩阵（对称正定）b=[12;-1;5];%常数向量%雅可比迭代法max_iter=1000;%最大迭代次数tol=1e-10;%收敛容差x6=zeros(size(b));%初始猜测值x_prev=x6;%上一次迭代值iter=0;%迭代计数器converged=false;%收敛标志%分解矩阵A=D+L+UD=diag(diag(A));%对角矩阵L=tril(A,-1);%下三角部分（不含对角线）U=triu(A,1);%上三角部分（不含对角线）%雅可比迭代公式：x^(k+1)=D^(-1)*(b-(L+U)*x^(k))whileiter<max_iter&&~convergedx6=D\(b-(L+U)*x_prev);error=norm(x6-x_prev)/norm(x6);%相对误差iferror<tolconverged=true;endx_prev=x6;iter=iter+1;end%验证解的正确性err6=norm(A*x6-b);%显示结果fprintf('线性方程组求解结果：\n');fprintf('----------------------------------------\n');fprintf('雅可比迭代：x=[%.6f,%.6f,%.6f],误差=%.10f\n',x6,err6);fprintf('雅可比迭代：迭代次数=%d,收敛=%f\n',iter,converged);fprintf('----------------------------------------\n');例2-9代码%线性方程组求解示例clearall;closeall;clc;%定义线性方程组Ax=bA=[4,-1,1;-1,4,-2;1,-2,4];%系数矩阵（对称正定）b=[12;-1;5];%常数向量%高斯-赛德尔迭代法max_iter=1000;%最大迭代次数tol=1e-10;%收敛容差x7=zeros(size(b));%初始猜测值x_prev=x7+inf;%上一次迭代值iter=0;%重置迭代计数器converged_gs=false;%高斯-赛德尔迭代收敛标志%高斯-赛德尔迭代公式：(D+L)x^(k+1)=b-Ux^(k)whileiter<max_iter&&~converged_gs%逐分量更新fori=1:length(x7)x7(i)=(b(i)-A(i,[1:i-1,i+1:end])*x7([1:i-1,i+1:end]))/A(i,i);enderror=norm(x7-x_prev)/norm(x7);%相对误差iferror<tolconverged_gs=true;endx_prev=x7;iter=iter+1;enditer_gs=iter;%保存高斯-赛德尔迭代次数%验证解的正确性err7=norm(A*x7-b);%显示结果fprintf('线性方程组求解结果：\n');fprintf('----------------------------------------\n');fprintf('高斯-赛德尔迭代：x=[%.6f,%.6f,%.6f],误差=%.10f\n',x7);fprintf('高斯-赛德尔迭代：迭代次数=%d,收敛=%f\n',iter_gs,converged_gs);fprintf('----------------------------------------\n');例2-10代码%非线性方程二分法求解示例clearall;closeall;clc;%定义要求解的非线性方程f=@(x)x.^3-2*x-5;%示例方程：f(x)=x³-2x-5%绘制函数图像，帮助确定根的大致位置x_plot=linspace(-3,3,1000);y_plot=f(x_plot);figure('Position',[100,100,800,600]);plot(x_plot,y_plot,'b-','LineWidth',2);holdon;plot(x_plot,zeros(size(x_plot)),'k--','LineWidth',1);%x轴title('函数f(x)=x³-2x-5的图像');xlabel('x');ylabel('f(x)');gridon;%从图像中观察，根位于区间[2,3]内a=2;%区间左端点b=3;%区间右端点tol=1e-6;%收敛容差max_iter=100;%最大迭代次数%检查区间端点的函数值符号是否相反iff(a)*f(b)>=0error('区间端点的函数值必须异号！');end%初始化迭代iter=0;c_prev=a;%上一次迭代的中点converged=false;%存储迭代过程iterations=zeros(max_iter,4);%存储[迭代次数,a,b,c]%二分法迭代fprintf('二分法迭代过程：\n');fprintf('迭代次数\t左端点(a)\t右端点(b)\t中点(c)\t\tf(c)\n');fprintf('------------------------------------------------------------------------\n');whileiter<max_iter&&~convergedc=(a+b)/2;%计算中点fc=f(c);%计算中点的函数值%存储迭代信息iterations(iter+1,:)=[iter+1,a,b,c];%输出当前迭代信息fprintf('%d\t\t%.8f\t%.8f\t%.8f\t%.8e\n',iter+1,a,b,c,fc);%检查收敛条件ifabs(fc)<tol||abs(c-c_prev)/abs(c)<tolconverged=true;end%更新区间iffc*f(a)<0b=c;%根在左半区间elsea=c;%根在右半区间endc_prev=c;iter=iter+1;end%输出结果fprintf('------------------------------------------------------------------------\n');fprintf('二分法求解结果：\n');fprintf('根的近似值：x=%.10f\n',c);fprintf('函数值：f(x)=%.10e\n',f(c));fprintf('迭代次数：%d\n',iter);fprintf('是否收敛：%f\n',converged);fprintf('区间长度：%.10e\n',b-a);%在图像上标记根的位置holdon;plot(c,0,'ro','MarkerSize',10,'MarkerFaceColor','r');legend('函数图像','x轴','根的位置');%绘制收敛过程figure('Position',[100,100,800,600]);iterations=iterations(1:iter,:);%截取有效数据subplot(2,1,1);semilogy(1:iter,abs(f(iterations(:,4))),'bo-','LineWidth',1.5);title('二分法收敛过程');xlabel('迭代次数');ylabel('|f(c)|(对数尺度)');gridon;subplot(2,1,2);semilogy(1:iter,(b-a)./2.^(1:iter),'r--','LineWidth',1.5);holdon;semilogy(1:iter,abs(iterations(:,4)-iterations(:,3)),'bo-','LineWidth',1.5);title('区间长度与误差估计');xlabel('迭代次数');ylabel('误差(对数尺度)');legend('理论误差','实际误差');gridon;axistight;例2-11代码%非线性方程求解方法对比clearall;closeall;clc;%定义要求解的非线性方程f=@(x)x.^3-2*x-5;%示例方程：f(x)=x³-2x-5df=@(x)3*x^2-2;%f(x)的导数%绘制函数图像，帮助确定根的大致位置x_plot=linspace(1,3,1000);y_plot=f(x_plot);figure('Position',[100,100,800,600]);plot(x_plot,y_plot,'b-','LineWidth',2);holdon;plot(x_plot,zeros(size(x_plot)),'k--','LineWidth',1);%x轴title('函数f(x)=x³-2x-5的图像');xlabel('x');ylabel('f(x)');gridon;%求解参数设置x0=2.5;%初始猜测值（用于迭代法）a=2;%二分法区间左端点b=3;%二分法区间右端点tol=1e-6;%收敛容差max_iter=100;%最大迭代次数%%方法1：二分法a_bisect=a;b_bisect=b;iter_bisect=0;c_prev=a_bisect;converged_bisect=false;whileiter_bisect<max_iter&&~converged_bisectc=(a_bisect+b_bisect)/2;fc=f(c);ifabs(fc)<tol||abs(c-c_prev)/abs(c)<tolconverged_bisect=true;endiffc*f(a_bisect)<0b_bisect=c;elsea_bisect=c;endc_prev=c;iter_bisect=iter_bisect+1;endx_bisect=c;err_bisect=abs(f(x_bisect));%%方法2：一般迭代法（不动点迭代）%将f(x)=0转换为x=g(x)的形式%例如：x=(2x+5)^(1/3)g=@(x)(2*x+5).^(1/3);x_iter=x0;x_prev=x_iter;iter_iter=0;converged_iter=false;whileiter_iter<max_iter&&~converged_iterx_iter=g(x_prev);error=abs(x_iter-x_prev);iferror<tol||abs(f(x_iter))<tolconverged_iter=true;endx_prev=x_iter;iter_iter=iter_iter+1;enderr_iter=abs(f(x_iter));%%方法3：牛顿迭代法x_newton=x0;x_prev=x_newton;iter_newton=0;converged_newton=false;whileiter_newton<max_iter&&~converged_newtonx_newton=x_prev-f(x_prev)/df(x_prev);error=abs(x_newton-x_prev);iferror<tol||abs(f(x_newton))<tolconverged_newton=true;endx_prev=x_newton;iter_newton=iter_newton+1;enderr_newton=abs(f(x_newton));%%结果对比fprintf('\n-------------------结果对比-------------------\n');fprintf('方法\t\t\t根的近似值\t\t迭代次数\t收敛\t\t误差\n');fprintf('----------------------------------------------------------------------------\n');fprintf('二分法\t\t\t%.10f\t\t%d\t\t%f\t\t%.2e\n',x_bisect,iter_bisect,converged_bisect,err_bisect);fprintf('一般迭代法\t\t%.10f\t\t%d\t\t%f\t\t%.2e\n',x_iter,iter_iter,converged_iter,err_iter);fprintf('牛顿迭代法\t\t%.10f\t\t%d\t\t%f\t\t%.2e\n',x_newton,iter_newton,converged_newton,err_newton);%%可视化三种方法的收敛过程figure('Position',[100,100,1200,800]);%绘制函数和根subplot(2,2,1);plot(x_plot,y_plot,'b-','LineWidth',2);holdon;plot(x_plot,zeros(size(x_plot)),'k--','LineWidth',1);plot(x_bisect,0,'ro','MarkerSize',10,'MarkerFaceColor','r','DisplayName','二分法');plot(x_iter,0,'go','MarkerSize',10,'MarkerFaceColor','g','DisplayName','一般迭代法');plot(x_newton,0,'mo','MarkerSize',10,'MarkerFaceColor','m','DisplayName','牛顿迭代法');title('函数图像与根的位置');xlabel('x');ylabel('f(x)');gridon;legend;%绘制收敛过程对比subplot(2,2,2);semilogy(1:iter_bisect,abs(f(x_bisect)-f(x_plot(1:iter_bisect))),'r-','LineWidth',1.5,'DisplayName','二分法');holdon;semilogy(1:iter_iter,abs(f(x_iter)-f(x_plot(1:iter_iter))),'g-','LineWidth',1.5,'DisplayName','一般迭代法');semilogy(1:iter_newton,abs(f(x_newton)-f(x_plot(1:iter_newton))),'m-','LineWidth',1.5,'DisplayName','牛顿迭代法');title('收敛过程对比');xlabel('迭代次数');ylabel('误差(对数尺度)');gridon;legend;%牛顿法迭代过程可视化subplot(2,2,3);plot(x_plot,y_plot,'b-','LineWidth',2);holdon;plot(x_plot,zeros(size(x_plot)),'k--','LineWidth',1);%绘制前几次牛顿迭代的切线x_steps=zeros(iter_newton,1);x_steps(1)=x0;fori=1:min(iter_newton,5)%只显示前5次迭代x_prev=x_steps(i);y_prev=f(x_prev);slope=df(x_prev);%切线方程：y=slope*(x-x_prev)+y_prevx_tangent=linspace(x_prev-0.5,x_prev+0.5,100);y_tangent=slope*(x_tangent-x_prev)+y_prev;plot(x_tangent,y_tangent,'--','LineWidth',1);plot(x_prev,y_prev,'o','MarkerSize',6,'MarkerFaceColor','r');%计算下一个迭代点ifi<iter_newtonx_next=x_prev-y_prev/slope;x_steps(i+1)=x_next;plot([x_prev,x_next],[y_prev,0],'r-','LineWidth',1);endendplot(x_newton,0,'mo','MarkerSize',10,'MarkerFaceColor','m');title('牛顿迭代法的几何解释');xlabel('x');ylabel('f(x)');gridon;%一般迭代法迭代过程可视化subplot(2,2,4);x_range=linspace(1.5,3,100);plot(x_range,x_range,'k-','LineWidth',1);%y=x直线holdon;plot(x_range,g(x_range),'b-','LineWidth',2);%y=g(x)曲线%绘制迭代路径x_steps=zeros(iter_iter,1);x_steps(1)=x0;y_steps=zeros(iter_iter,1);y_steps(1)=g(x0);fori=1:min(iter_iter,10)%只显示前10次迭代x_prev=x_steps(i);y_prev=y_steps(i);%绘制从(x_prev,y_prev)到(x_next,y_prev)的水平线ifi<iter_iterx_next=y_prev;y_next=g(x_next);x_steps(i+1)=x_next;y_steps(i+1)=y_next;%绘制水平线和垂直线plot([x_prev,x_next],[y_prev,y_prev],'r-','LineWidth',1);plot([x_next,x_next],[y_prev,y_next],'r-','LineWidth',1);endend%标记初始点和不动点plot(x0,g(x0),'ro','MarkerSize',6,'MarkerFaceColor','r');plot(x_iter,g(x_iter),'mo','MarkerSize',10,'MarkerFaceColor','m');title('一般迭代法的几何解释(y=g(x))');xlabel('x');ylabel('g(x)');gridon;legend('y=x','y=g(x)');axisequal;有限元代码11：主程序clcclear;node=[1000;2100;3110;4010;];%节点信息，第一列为节点编号，2~4列分别为x,y,z方向坐标ele=[1123;2341;];%单元信息，第一列为单元编号，后面各列为单元上的节点号码num_ele=size(ele,1);%单元数%---物理参数------------------E=2.1e11;%弹性模量MPat=0.01;%单元厚度mmiu=1/3;%泊松比q=1e6;%均布载荷N/m%-------------------------------n_ele=length(ele(:,1));%单元数%组装总体刚度矩阵dof=length(node(:,1))*2;%自由度数f=ones(dof,1)*1e8;%结构整体外载荷矩阵，整体坐标系下f_loc=zeros(6,1);%单元外载荷矩阵，局部坐标系下u=ones(dof,1)*1e6;%位移矩阵K=zeros(dof);%总体刚度矩阵stress=zeros(n_ele,1);%单元应力矩阵fori=1:n_elek_ele=TriangleElementStiffness(E,miu,t,node(ele(i,2:4),2:4));K=assemTriangle(K,k_ele,ele(i,2),ele(i,3),ele(i,4));end%力边界条件f(6)=q/2;%3节点垂向力f(8)=q/2;%4节点垂向力f(3)=0;f(5)=0;%位移边界条件u(1)=0;u(2)=0;u(4)=0;u(7)=0;%求解未知自由度index=[];%未知自由度的索引p=[];%未知自由度对应的节点力矩阵fori=1:dofifu(i)~=0index=[index,i];p=[p;f(i)];endendu(index)=K(index,index)\p;%高斯消去f=K*u;%单元应力stress=zeros(num_ele,3);x1=node(:,2)+u(1:2:end);%变形后的位移y1=node(:,3)+u(2:2:end);%变形后的位移figure;fori=1:n_eleu1=[u(2*ele(i,2)-1);u(2*ele(i,2));u(2*ele(i,3)-1);u(2*ele(i,3));u(2*ele(i,4)-1);u(2*ele(i,4))];stress(i,:)=TriangleElementStress(E,miu,node(ele(i,2:4),2:3),u1,1)';%单元应力计算patch(node(ele(i,2:4),2),node(ele(i,2:4),3),stress(i,1));endholdon;figure;fori=1:n_elepatch(node(ele(i,2:4),2),node(ele(i,2:4),3),'w','FaceColor','none','LineStyle','-','EdgeColor','b');%变形前holdon;patch(x1(ele(i,2:4)),y1(ele(i,2:4)),'w','FaceColor','none','EdgeColor','r');%变形后end2：单元刚度子函数functionk_ele=TriangleElementStiffness(E,miu,t,node_ele)%TriangleElementStiffnessThisfunctionreturnstheelement%stiffnessmatrixforaTriangle(CST)%elementwithmodulusofelasticityE,%Poission'sratiomiu,constantthicknesst,%node_elethenodecoordinateofelement.%Thesizeoftheelementstiffness%matrixis6x6.%---------nodecoordinate------x1=node_ele(1,1);y1=node_ele(1,2);x2=node_ele(2,1);y2=node_ele(2,2);x3=node_ele(3,1);y3=node_ele(3,2);%-------------------------------A=(x1*(y2-y3)+x2*(y3-y1)+x3*(y1-y2))/2;%单元面积%A=abs(A);a1=x2*y3-y2*x3;a2=y1*x3-x1*y3;a3=x1*y2-y1*x2;b1=y2-y3;b2=y3-y1;b3=y1-y2;c1=x3-x2;c2=x1-x3;c3=x2-x1;B=1/2/A*[b10b20b30;0c10c20c3;c1b1c2b2c3b3];D=E/(1-miu^2)*[1miu0;miu10;00(1-miu)/2];%strss/strainmatrixforplanestress%D=E/(1-miu^2)*[1miu0;%miu10;%00(1-miu)/2];%strss/strainmatrixforplanestraink_ele=t*A*B'*D*B;%的单元刚度矩阵3：单元应力子函数functionstr=TriangleElementStress(E,miu,node_ele,u1,p)%TriangleElementStressThisfunctionreturnstheelement%stressmatrixforaTriangle%elementwithmodulusofelasticityE,%Poission'sratiomiu,%node_elethenodecoordinateofelement，planestressorplanestrainoption.%---------nodecoordinate------x1=node_ele(1,1);y1=node_ele(1,2);x2=node_ele(2,1);y2=node_ele(2,2);x3=node_ele(3,1);y3=node_ele(3,2);%-------------------------------A=(x1*(y2-y3)+x2*(y3-y1)+x3*(y1-y2))/2;%单元面积a1=x2*y3-y2*x3;a2=y1*x3-x1*y3;a3=x1*y2-y1*x2;b1=y2-y3;b2=y3-y1;b3=y1-y2;c1=x3-x2;c2=x1-x3;c3=x2-x1;B=1/2/A*[b10b20b30;0c10c20c3;c1b1c2b2c3b3];ifp==1D=E/(1-miu^2)*[1miu0;miu10;00(1-miu)/2];%strss/strainmatrixforplanestresselseifp==2D=E/(1+miu)/(1-2*miu)*[1-miumiu0;miu1-miu0;00(1-2*miu)/2];%strss/strainmatrixforplanestrainendstr=D*B*u1;%单元应力4：整刚组装子函数functionk_t=assemTriangle(k_t,k_ele,node1,node2,node3)%assemTriangleThisfunctionassemblestheelementstiffness%matrixkoftheplaneTriangleelementwithnodes%iandjintotheglobalstiffnessmatrixK.%Thisfunctionreturnstheglobalstiffness%matrixKaftertheelementstiffnessmatrix%kisassembled.d(1:2)=2*node1-1:2*node1;d(3:4)=2*node2-1:2*node2;d(5:6)=2*node3-1:2*node3;forii=1:6forjj=1:6k_t(d(ii),d(jj))=k_t(d(ii),d(jj))+k_ele(ii,jj);endend有限元代码2%有限元程序，非子