科学计算与仿真应用（基于MWORKS）课件 第4、5章 Syslab绘图功能教学、数值计算方法

第4章Syslab绘图功能教学本章目标了解Syslab中绘图功能的整体框架掌握

二维图形的绘制方法与常用函数（plot等）熟悉

三维图形绘制的基本流程（plot3、mesh等）学会

图形属性设置方法（坐标轴、网格、标题、图例等）能够根据实际问题

选择合适的绘图函数并完成可视化表达知识点总览二维绘图函数：掌握plot及其多种调用形式，实现基本曲线绘制。图形属性设置：学会设置坐标轴、网格线、标题、图例等图形元素。特殊二维图形：了解条形图、直方图、极坐标图等常见可视化方法。简易绘图函数：掌握ezplot、ezpolar等快速绘图工具。三维绘图基础：理解plot3、patch等三维图形绘制方法。三维曲面与网格：学会mesh、meshc、meshz等曲面可视化技术。plot函数:二维绘图的核心plot函数是Syslab中用于绘制二维图形的基础函数,支持多种数据类型和丰富的属性设置基本语法格式plot(x)实向量→以下标为横坐标;实矩阵→以列下标为横坐标;复数矩阵→实部为横坐标、虚部为纵坐标plot(x,y)同维数组→以元素为横纵坐标;一方为数组、另一方为矩阵→绘制多条曲线plot(x1,y1,x2,y2)在同一窗口绘制多条曲线线色、点符号、线型线色b蓝g绿r红c青m洋红y黄k黑w白点符号.点o圆x叉+加号*星号D钻石v倒三角线型-实线:点线-.点划线--虚线plot函数:二维绘图的核心其他常用属性linewidth线宽设置markersize点符号大小markerfacecolor点填充颜色markeredgecolor点边框颜色代码示例usingTyPlotusingTyBasex

=

0:

2*pi/30

:

2*piy

=

sin.(x)plot(x,

y,

"ro-.",

linewidth=1,

markersize=10,

markerfacecolor="g",

markeredgecolor="y")绘制红色虚线,带圆形标记,线宽为1,标记大小为10,填充绿色,边框黄色图形参数设置通过图形参数设置,可以精确控制坐标轴、标题、分格线等元素,提升图表的专业性和可读性分格线控制grid("on")

#

显示分格线grid("off")

#

隐藏分格线grid

#

切换状态分格线的疏密程度由坐标刻度决定,如需改变需先定义坐标刻度标题与标签title("文字")添加图形标题,支持fontname和fontsize属性xlabel/ylabel设置坐标轴标签legend("字符串")添加图例说明text(X,Y,"字符串")在指定位置添加文本坐标轴设置axis([xminxmaxyminymax])设置x轴和y轴显示范围axis("auto")返回坐标刻度默认值axis("equal")各坐标轴等刻度显示axis("tight")坐标范围与数据区间一致axis("off")/axis("on")关闭/打开坐标系图形参数设置子窗口设置subplot(m,n,p)将窗口拆分为m×n个子窗口,第p个为当前窗口子窗口顺序:先按行从左到右特殊符号支持希腊字母\alphaα\betaβ\gammaγ\thetaθ\piπ\omegaω\sigmaσ\phiφ数学符号\infty∞\int∫\partial∂\rightarrow→图例与注释legend("说明")：添加图例text(x,y,"字符串")：在图中插入说明文字作用：区分多条曲线、突出关键信息特殊二维图形Syslab提供多种特殊二维图形函数,满足不同数据可视化需求条形图与直方图bar(x,y)竖直条形图,x为横坐标,y为向量或矩阵barh(x,y)水平条形图hist(y,n)直方图,n为区间数量特殊线图stairs(x,y)阶梯图形stem(x,y)茎状图plt_fill(x,y,c)填充图,c指定填充颜色极坐标与矢量图polarplot(theta,rho)极坐标图,theta为极角,rho为径向长度feather(u,v)矢量图,u为x轴速度,v为y轴速度其他实用函数fplot快速绘图comet彗星轨迹pie饼图scatter散点图特殊二维图形Easy绘图函数简化版绘图函数,无需精确定义数据点ezplot(fun,[xmin,xmax])绘制函数在指定域内的图形ezplot(funx,funy,[tmin,tmax])绘制参数方程曲线ezpolar(fun,[a,b])绘制极坐标曲线Easy函数适合快速查看函数图形,不需要预先计算数据点特殊二维图形的核心作用：让“数据结构”比“数值本身”更直观连续变化

→

曲线图分类对比

→

条形图分布特征

→

直方图方向信息

→

极坐标或矢量图通常来说，选择图形的原则如下：plot3与patch:三维绘图基础三维绘图从plot3函数开始,通过patch函数可以创建填充多边形plot3函数plot3(x,

y,

z,

s)x、y、z:三维坐标数据,尺寸必须相等s:设置线型、颜色、数据点标记的字符串支持在同一窗口绘制多条三维曲线(使用hold("on"))主要功能绘制空间轨迹展示参数曲线支持多条曲线叠加代码示例-螺旋曲线usingTyPlotusingTyBaset

=

LinRange(0,12*pi,100)plot3(t.*sin.(t),

t.*cos.(t),

t,

"b-o",

markerfacecolor="g",

linewidth=2)xlabel("sin(t)")ylabel("cos(t)")zlabel("t")grid("on")plot3与patch:三维绘图基础patch函数patch(x,

y,

z,

c)用途:建立补片对象(多边形)c:指定填充颜色多边形未封闭时会自动封闭主要功能构建几何模型绘制立体结构表示空间区域代码示例-绘制有填充面的立方体usingTyPlotusingTyBasefigure()#打开一个新图形窗口.#定义x1、y1、z1、x2、y2、z2、x3、y3、z3，表示三个填充面的顶点坐标.x1=[0,1,1,0];y1=[0,0,0,0];z1=[0,0,1,1];x2=[0,1,1,0];y2=[0,0,1,1];z2=[1,1,1,1];x3=[0,0,0,0];y3=[0,0,1,1];z3=[0,1,1,0];#填充由三个数组x、y、z定义的区域.patch(x1,y1,z1,"y");hold("on")patch(x2,y2,z2,"b");patch(x3,y3,z3,"g");plot3与patch:三维绘图基础应用场景plot3应用三维曲线、螺旋线、空间轨迹、参数方程曲线patch应用立方体、多面体、填充曲面、三维建模基础主要区别plot3表达“线结构”（轨迹、路径）patch表达“面结构”（区域、几何体）三维网格图与曲面图mesh和surf函数是三维可视化中最常用的函数,分别用于绘制网格图和曲面图peaks函数生成测试数据,常用于mesh、surf等函数x,y,z

=

peaks(n)生成n×n矩阵,默认n=49表达式包含高斯函数的组合shading阴影模式shadingflat平面阴影模式shadingfaceted面元阴影模式(默认)mesh函数(网格图)mesh(x,y,z,c)基于x、y、z矩阵绘制网格图hidden("on")/hidden("off"):控制是否显示被遮挡部分meshc:带等高线的网格图meshz:带帘幕的镂空网格图surf函数(曲面图)surf(x,y,z,c)绘制颜色填充的网格图与mesh的区别:surf填充面元surfc:带等高线的曲面图surfnorm:带法线的曲面图三维网格图与曲面图颜色设置colormap(map)设置色图,常用选项:jethsvhotcoolspringsummerautumnwintercolorbar

#

显示颜色条代码示例-绘制峰面usingTyPlotx,y,z=peaks(30);figure()surf(x,y,z);colormap("summer")title("峰面")等值线与矢量图contour函数用于绘制等值线,quiver函数用于绘制矢量场,是科学计算可视化的重要工具contour函数(等值线)contour(z)绘制z的等值线图contour(x,y,z)指定x、y轴坐标contour(z,n)n设置等值线条数contour(z,v)v向量定义等值线数值c,h=contour(...)返回等值线矩阵c和句柄hcontour3与contourfcontour3(x,y,z)绘制三维等值线图contour3(x,y,z,n)设置等值线条数contourf(z,n)绘制填充的等值线图clabel(c)为等值线添加数值标注等值线与矢量图quiver函数(矢量图)quiver(x,

y,

u,

v)绘制由小箭头构成的矢量图起始点:(x,y)对应的格点箭头长度:由(u,v)值确定代码示例x,y,z

=

peaks(30)vx,vy

=

gradient(z,2,2)contour(x,y,z,10)hold("on")quiver(x,y,vx,vy)axis("image")结合contour和quiver,可在等值线图上叠加矢量场综合应用实例实例1:带洞的峰面利用NaN(NotaNumber)在曲面上创建"洞"usingTyBase#加载TyBase函数库.#借助peaks函数定义变量x、y、z.x,y,z=peaks(50);#借助NaN（Notanumber）在峰面上定义一个洞.z[20:30,15:25].=NaN;figure()#打开一个新图形窗口.#绘制带洞的峰面.a=surf(x,y,z)#设置colormap属性为winter.colormap(a,"winter");#设置shading属性为faceted.shading(a,"faceted")将特定区域的z值设为NaN,即可在曲面上创建"洞"综合应用实例实例2:二重透视球面使用hold("on")在同一窗口绘制多个球面,通过alpha设置透明度实现透视效果#借助sphere函数定义矩阵x、y、z.x,y,z=sphere(30;fig=false)#定义x1、y1、z1和x2、y2、z2两套值.x1=15.*x;y1=15.*y;z1=15.*z;x2=20.*x;y2=20.*y;z2=20.*z;figure()#打开一个图形窗口.#绘制最小的球面.a=mesh(10.*x,10.*y,10.*z);#保留在当前图形窗口.hold("on")alpha(0)设置透明,可看到内部结构#绘制第二个球面.b=mesh(x1,y1,z1);#设置colormap属性为autumn.colormap(b,"autumn")c=mesh(x2,y2,z2);#绘制最大的球面.alpha(0)#设置透明度属性为0.#设置colormap属性为summer.colormap(c,"summer")hold("off")#释放当前图形窗口.axis("square")#设置axis属性为square.axis("off")#隐藏坐标系.综合应用实例实例3:缺失网格透视图结合NaN和colormap创建特殊视觉效果不同球面使用不同色图autumn、summer等色图部分区域设为NaN图形属性鼠标操作点击工具栏图标右击图形打开属性菜单设置面颜色、边颜色、线型、线宽"属性检查器"提供更多设置#借助sphere函数定义矩阵x、y、z.x,y,z=sphere(30;fig=false)#定义x1、y1、z1和x2、y2、z2两套值.x1=1.5.*x;y1=1.5.*y;z1=1.5.*z;x2=2.*x;y2=2.*y;z2=2.*z;z1[1:31,3:13].=NaN;z2[1:31,3:13].=NaN;#绘制最小的球面.a=surf(x,y,z);shading(a,"flat");#保留在当前图形窗口.hold("on")#绘制第二个球面.b=mesh(x1,y1,z1);#设置colormap属性为autumn.colormap(b,"autumn")#绘制最大的球面.c=mesh(x2,y2,z2);#设置colormap属性为summer.colormap(c,"summer")hold("off")#释放当前图形窗口.axis("square")#设置axis属性为square.axis("off")#隐藏坐标系.本章小结本章系统介绍了Syslab绘图相关函数及其应用方法，重点围绕二维与三维数据可视化展开。（1）二维绘图是数据可视化的基础，本章以plot函数为核心，介绍了其多种调用形式及参数设置方法，包括向量、矩阵和复数数据的绘制规则，并说明了如何在同一坐标系中绘制多条曲线，为后续图形表达奠定基础。（2）为了使图形表达更加规范和清晰，本章详细介绍了图形参数设置方法，包括网格线控制、坐标轴范围设置、标题与坐标轴标签、图例说明以及文本注释等内容，使绘图结果不仅能够展示数据，还能够准确传达物理意义。（3）在基本曲线绘制的基础上，本章进一步介绍了多种特殊二维图形函数，如条形图、直方图、茎状图、阶梯图、极坐标图等，这些图形能够更直观地反映数据分布、离散特征及方向信息，从而丰富数据表达方式。（4）在三维可视化部分，本章介绍了plot3与patch等基础函数，实现空间曲线与几何结构的绘制，并结合网格图与曲面图函数说明三维数据的展示方法，使读者能够从平面数据分析过渡到空间数据可视化。通过本章学习，读者可以建立从二维曲线绘制到三维图形表达的完整绘图框架，掌握利用可视化手段分析和展示数据的基本方法，为后续建模分析与工程应用提供直观有效的工具。第5章数值计算方法本章目标理解并掌握在Syslab中Julia语言下的线性方程组的直接法与迭代法，以及矩阵特征值和特征向量的常用数值方法（如幂法、反幂法、Jacobi方法和QR方法）。掌握syslab中多项式插值与最小二乘曲线拟合的基本原理，能够应用这些方法进行数据拟合与函数逼近。理解数值积分与数值微分的基本方法，掌握常微分方程的Euler方法和Runge-Kutta方法，并了解相关函数的使用。同时，掌握使用julia语言对非线性方程求根的常用数值方法，如二分法、Newton迭代法和弦截法。知识点总览线性方程组的数值解法：掌握直接法（如LU分解、Cholesky分解、矩阵除法）和迭代法（Jacobi迭代、Gauss-Seidel迭代、SOR迭代）的基本原理与实现，能够针对不同类型（稠密/稀疏、正定/非正定）的线性方程组选择合适的数值解法。多项式插值与数据拟合技术：理解Lagrange插值、Hermite插值、样条插值等方法，掌握最小二乘曲线拟合的原理与实现，能够对离散数据进行函数逼近与预测。数值积分与数值微分的计算方法：掌握自适应Simpson积分、Romberg算法等数值积分方法，以及基于插值多项式的数值微分公式（如三点公式、差分法），能够对复杂函数或离散数据进行积分与导数的近似计算。矩阵特征值与特征向量的数值求解：理解幂法、反幂法、Jacobi方法、QR方法等经典算法的基本思想与适用场景，能够编写程序求解矩阵的主特征值、全部特征值及对应的特征向量。常微分方程与非线性方程的数值解法：掌握常微分方程初值问题的Euler方法、改进Euler方法、Runge-Kutta方法，以及非线性方程求根的二分法、Newton迭代法、弦截法，能够分析算法的精度、收敛性与稳定性。1.线性方程组的数值解法-直接法#将解向量中的元素转换为分数形式.x_fraction=map(x->Rational(x),x)直接法是经过有限步的算术运算求线性方程组精确解的方法（假设计算过程中没有舍入误差）。但由于在实际计算过程中存在舍入误差，因此这种方法只能求得线性方程组的近似解。示例：求解线性方程组运算规则：对于一个n阶线性方程组Ax=b，在Syslab中，只需要一个反斜杠运算符“/”或斜杠运算符“\”就可求解。#对矩阵A和b赋值.A=[5-41;-46-4;1-46];b=[2;-1;-1]x=A\b1.线性方程组的数值解法-迭代法#使用反斜杠运算符求解x.result=Jacobimethod(A,b,x0,Maxit);示例：使用Jacobi迭代法求解如下线性方程组，要求Jacobi迭代法：迭代的收敛性只和系数矩阵有关，和初值的选择没有关系，初值可任意选择。时停止：注意：需要先编辑Jacobi函数，再代入数值1.线性方程组的数值解法-迭代法#SOR更新公式.x1[i]=(1-Omega)*x0[i]+(Omega/A[i,i])*(b[i]-_sum1-_sum2);示例：求解如下线性方程组，要求SOR（SuccessiveOver-Relaxation，逐次超松弛）迭代法：是Gauss-Seidel迭代的一种加速方法，是解大型稀疏矩阵线性方程组的有效方法之一，其计算公式简单，但需要较好的加速因子（最佳松弛因子）。时停止：注意：SOR迭代法中最佳松弛因子ω的选择很困难，数值分析已经证明SOR迭代收敛的必要条件是ω在0和2之间，并且在系数矩阵是三角矩阵的情况下，给出了一个ω的计算公式。2.多项式插值与最小二乘曲线拟合-多项式插值Lagrange插值：已知个插值点及其对应的函数值若采用Lagrange插值方法，则对插值区间内任意x可以求得示例：已知数据如表5.1所示。试使用Lagrange插值方法求时，函数的近似值。

表5.1例5.6数据0.561600.562800.564010.565210.827410.826590.825770.814952.多项式插值与最小二乘曲线拟合-多项式插值functionlagrange(x,y,x0)n=length(x)m=length(x0)y0=zeros(m)k=0;i=0;j=0;forkin1:mz=x0[k];s=0.0;foriin1:np=1.0;forjin1:nifj!=ip*=(z-x[j])/(x[i]-x[j]);endends+=p*y[i];endy0[k]=s;endreturny0end使用Julia语言编写lagrange函数，之后赋值进行计算2.多项式插值与最小二乘曲线拟合-多项式插值代入数据后得到运行结果：pchip函数和spline函数的插值结果如图5.2所示。图5.22.多项式插值与最小二乘曲线拟合-最小二乘曲线拟合最小二乘曲线拟合：具体做法是，对于给定的一组数据functionlinearfit(x,y,m)nx=length(x)ny=length(y)ifnx!=nyerror("Theinputdataiswrong!");endA=zeros(m+1,m+1);b=zeros(m+1);foriin1:m+1forjin1:m+1A[i,j]=sum(x.^(i+j-2));endb[i]=sum(x.^(i-1).*y);enda=A\b;yy=polyval(reverse(a),x)，要求在函数类中找到一个函数，使误差的平方和最小，即温度/℃100110120130140150160170180190产品率/％45515461667074788589示例：已知某种化学反应过程中温度与产品率关系的数据如表5.2所示。请用最小二乘曲线拟合得出函数关系。2.多项式插值与最小二乘曲线拟合-最小二乘曲线拟合小二乘曲线拟合结果如图5.3所示图5.33.数值积分与数值微分-数值积分（1）quad函数。quad函数是Syslab自带的函数，其功能是使用自适应Simpson公式求给定函数的定积分。#q=quad(fun,a,b);q=quad(fun,a,b;Name=Value);语法规则：根据节点的选择方法，数值分析中的数值积分方法有Newton-Cotes公式法、Romberg算法、Gauss公式法等。示例：求积分注意：在运行脚本前需要安装QuadGK函数库，其语法格式如下：#julia>importPkg;julia>Pkg.add("QuadGK")usingTyMathusingQuadGKf(x)=sqrt(4-sin(x)^2)result,error=quadgk(f,0,pi/6)println("积分结果:",result)println("估计误差:",error)3.数值积分与数值微分-数值积分（2）Romberg算法和函数。usingTyMathfunctionromberg_integration(f,a,b,n=6)R=zeros(n,n);h=(b-a)./(2.^(0:n-1));R[1,1]=(b-a)*(f(a)+f(b))/2;foriin2:nsum_mid=0.0;step=2^(n-i);forjin1:2^(i-2)x=a+(2*j-1)*step*h[i];sum_mid+=f(x);endR[i,1]=R[i-1,1]/2+sum_mid*h[i];forkin2:iR[i,k]=R[i,k-1]+(R[i,k-1]-R[i-1,k-1])/(4^(k-1)-1);endendintegral=R[n,n];error=abs(R[n,n]-R[n,n-1]);returnintegral,errorendclear()clc()f(x)=sqrt(4-sin(x)^2)result,error=romberg_integration(f,0,pi/6,6)println("Lundberg积分结果:",result)println("估计误差:",error)3.数值积分与数值微分-数值微分本节讨论数值微分的方法。插值型求导公式用插值多项式的导数近似表示，即#functionnpointdiff(x,y,x0)A=copy(y);N=length(x);forjin2:NforkinN:-1:j#Calculatethej-thorderdifferenceatpointk.A[k]=(A[k]-A[k-1])/(x[k]-x[k-j+1]);endendX=x[1];df=A[2];prod=1;n1=length(A)–1;forkin2:n1prod*=(x0-x[k])df+=prod*A[k+1];endreturnA,dfend基于

个点的差分求解程序如下。示例：已知函数的数据如表5.4所示，计算处的一阶导数。x2.52.62.72.82.9y12.18513.463714.879716.444618.1741usingTyMathusingTyBasex=[2.5,2.6,2.7,2.8,2.9];y=[12.185,13.4637,14.8797,16.4446,18.1741];A,df=npointdiff(x,y,2.7);println("A是:",A)returnAprintln("df是:",df)returndf4.矩阵的特征值和特征向量

Syslab中求解特征值和特征向量的函数如下：eig函数：求矩阵的特征值和特征向量。eigvals函数：求矩阵的特征值。eigvecs函数：求矩阵的特征向量。svd函数：计算矩阵的奇异值分解。schur函数：Schur分解。示例：求矩阵的特征值和特征向量。clear()clc()usingTyMathusingTyBaseA=[2-10;-12-1;0-12];V,D=eig(A);println("矩阵A的特征向量为：",V)returnVprintln("矩阵A的特征值为：",D)returnD4.矩阵的特征值和特征向量-幂法和反幂法functionpowermethod(A,x0,Eps=1e-6,maxtimes=500)Lambda=0;cnt=0;err=1;state=1;while(cnt<=maxtimes)&&(state==1)Y=A*x0;m,j=findmax(abs.(Y));c1=mdc=abs(Lambda-c1);Y=(1/c1)*Y;dv=norm(x0-Y);err=max(dc,dv);x0=Y;Lambda=c1;state=0;幂法是用来求矩阵的模最大的特征值及对应的特征向量的迭代法。迭代公式如下：使用幂法计算矩阵的模最大的特征值及对应的特征向量的自编函数如下：iferr>Epsstate=1endcnt+=1end#返回结果.returnLambda,x0end示例：使用幂法求矩阵的模最大的特征值及对应的特征向量。clear()clc()usingTyMathusingTyBaseA=[011-5;-217-7;-426-10];x0=[1;1;1]；Lambda,V=powermethod(A,x0,1e-6,500)println("矩阵A的模最大的特征值为：",Lambda)returnLambdaprintln("矩阵A的模最大的特征值对应的特征向量为：",V)returnV4.矩阵的特征值和特征向量-幂法和反幂法反幂法是用来求矩阵的模最小的特征值及对应的特征向量的迭代法。用替代执行上述的幂法就是所谓的反幂法。使用反幂法计算矩阵的模最小的特征值及对应的特征向量的自编函数如下：functioninvpowermethod(A,x0,alpha,Eps,maxtimes)

n,_=size(A);A=A-alpha*I;Lambda=0;cnt=0;err=1;state=1;while(cnt<=maxtimes)&&(state==1)Y=A\x0;m,j=findmax(abs.(Y));c1=m;dc=abs(Lambda-c1);Y=(1/c1)*Y;dv=norm(x0-Y);err=max(dc,dv);x0=Y;

Lambda=c1;state=0;iferr>Epsstate=1;endcnt+=1;endLambda=alpha+1/c1;V=x0;returnLambda,Vend

4.矩阵的特征值和特征向量-Jacobi方法Jacobi方法是用来计算实对称矩阵的全部特征值和特征向量的一种迭代法。使用Jacobi方法计算矩阵的特征值和特征向量的自编函数如下：functionjacobimethod(A,tol=1e-8,maxiter=1000)n=size(A,1);V=Matrix{Float64}(I,n,n);iters=0;whiletruemaxval=0.0;p,q=0,0;foriin1:n-1forjini+1:nifabs(A[i,j])>maxvalmaxval=abs(A[i,j]);p,q=i,j;endendendifmaxval<tol||iters>=maxiterbreakendθ=0.5*atan(2*A[p,q],A[q,q]-A[p,p]);R=Matrix{Float64}(I,n,n);c,s=cos(θ),sin(θ);R[p,p]=c;R[q,q]=c;R[p,q]=s;R[q,p]=-s;A=R'*A*R;V*=R;iters+=1;endidx=sortperm(diag(A));eigenvalues=diag(A)[idx];eigenvectors=V[:,idx];println("矩阵的特征值为:",eigenvalues)println("矩阵的特征向量为:",eigenvectors)end4.矩阵的特征值和特征向量-QR方法QR方法是目前求一般方阵的全部特征值最有效且被广泛应用的方法之一。通过一系列正交相似变换得到方阵序列usingLinearAlgebraXfunctionqrmethod(A::Matrix{<:Number};max_iter=1000,tol=1e-10)Ak=copy(A)n=size(Ak,1);foriterin1:max_iterQ,R=qr(Ak);Ak=R*Q;sub_diag_norm=norm([Ak[i+1,i]foriin1:n-1]);ifsub_diag_norm<tolbreakendendreturndiag(Ak)end{Ak}，在一定条件下，可以证明当，Ak时本质收敛于一个块对角矩阵，其对角子块是1阶或2阶方阵，因此容易求出它们的特征值。使用QR方法计算矩阵的特征值的自编函数如下。5.常微分方程的数值解法-Euler方法对于简单的一阶方程的初值问题，即functioneulermethod(f,x0,y0,h,n_steps)x_values=zeros(n_steps+1);y_values=zeros(n_steps+1);x_values[1]=x0;y_values[1]=y0;foriin1:n_stepsx=x_values[i];y=y_values[i];y_next=y+h*f(x,y);x_next=x+h;x_values[i+1]=x_next;y_values[i+1]=y_next;endreturnx_values,y_valuesend，由Euler方法可得使用Euler方法求解一阶方程的初值问题的自编函数如下。5.常微分方程的数值解法-Euler方法示例：使用Euler方法求解：并与精确解进行比较。根据自编程序代入数据后得到运行结果如图5.4所示：图5.45.常微分方程的数值解法-Runge-Kutta方法为了提高常微分方程的数值解法的精度，可间接使用Taylor展开的思想构造高阶的数值解法，四阶Runge-Kutta方法的经典公式如下：5.常微分方程的数值解法-Runge-Kutta方法使用Runge-Kutta方法求解常微分方程的自编函数如下。functionrungekuttamethod(f,a,b,y0,h)x=a:h:b;n=length(x);y=zeros(n);y[1]=y0foriin1:n-1k1=f(x[i],y[i]);k2=f(x[i]+h/2,y[i]+h/2*k1);k3=f(x[i]+h/2,y[i]+h/2*k2);k4=f(x[i+1],y[i]+h*k3);y[i+1]=y[i]+(k1+2*k2+2*k3+k4)*h/6;endreturnx,yendfunctionmyfun(x,y)return-1.0/xend5.常微分方程的数值解法表5.5例5.17的三种数值解法的对比i精确值Euler方法的误差改进的Euler方法的误差Runge-Kutta方法的误差01-100011.1-0.90910.0091-0.0009-0.000002405721.2-0.83330.0134-0.0013-0.000003416631.3-0.76920.0152-0.0015-0.000003723241.4-0.71430.0156-0.0015-0.000003676951.5-0.66670.0153-0.0015-0.00000346095.常微分方程的数值解法-Syslab中的相关函数函数算法适用条件精度特点ode454阶、5阶Runge-Kutta方法非刚性方程中最常用的单步算法。ode113可变阶算法非刚性方程低—高多步算法，在相同的精度下，它可能比ode45函数或ode23函数运算速度快一些，但不适用于不连续的系统。ode15s可变阶算法刚性方程低—中多步算法，方程刚性，用于ode45函数运算速度很慢时。ode232阶、3阶Runge-Kutta方法非刚性方程低有时对于轻度的刚性方程来说，它比ode45函数更有效。在相同的精度下，它需要比ode45函数更小的步长。ode23s基于改进的Rosenbrock公式刚性方程低单步算法，用于不适合使用ode15s函数的刚性方程。ode23t自由内插实现的梯形规则一般刚性低给出的解无数值衰减。ode23tbTR-BDF2算法刚性方程低当允许误差范围较大时，它比ode15s函数好。表5.6常微分方程的几种数值求解函数的对比6.非线性方程求根的数值计算方法-二分法

重复上述步骤，如果通过有限步求得根，则计算中止，否则每次有根区间长度减半，并得到一个闭区间：使用二分法求解非线性方程的自编函数如下。functionbisectionmethod(f,a,b;tol=1e-6,max_iter=100)iff(a)==0returna,0,trueendiff(b)==0returnb,0,trueendiff(a)*f(b)>0error("区间[a,b]两端点函数值必须异号")endconverged=false;iter=0;c=0.0;foriterin1:max_iterc=(a+b)/2;fc=f(c);ifabs(fc)<tol||(b-a)/2<tolconverged=true;break;endiffc*f(a)<0b=c;elsea=c;endendreturnc,iter,convergedend6.非线性方程求根的数值计算方法-Newton迭代法解非线性方程的Newton迭代法是一种将非线性函数线性化