实验七-常微分方程解法.ppt

1,数值实验六,数值积分与数值微分,2,实验报告要求,实验报告七必须在12月22日前上交算法使用编程语言不限；电子版试验报告与程序代码发送至邮箱：文件名称：实验七 08119000 张三,3,Matlab基础（实验一）,基本命令基本数学运算符号解函数文件格式函数调用文件保存位置,4,基本命令,1.help 帮助命令 help format%查看format命令的帮助；2.clear 清除变量命令 clear x%清除变量x(无参数时清理工作空间）3.清理命令窗口 clc%清屏4.format 输出格式 format long%15位有效数字格式输出 format short%4-7位有效数字格式输出 format compact%紧凑格式(输出不加空行）,5,基本数学运算,乘法：*除法：/乘方：根号：sqrt(x)正弦：sin(x)余弦：cos(x)自然对数：log(x)自然指数：exp(x),6,1.vpa 控制运算精度vpa(表达式,运算精度)vpa(pi,20)%显示至100位有效数字2.solve 方程的符号解solve(表达式,x),符号解,7,符号解,例 解方程解 在MATLAB工作窗口输入命令:y=solve(x3-sin(x)-12*x+1=0,x)vpa(y,5）y=roots(2,0,-1,-1),8,函数文件格式,函数文件由function语句引导，其格式为：function 输出形参表=函数名(输入形参表)%注释说明部分函数体：例：function k,xk,yk,p=jhnewtonqx(x0,ddmax)%牛顿切线法求非线性方程的根注：其中函数名的命名规则与变量名相同。输入形参为函数的输入参数，输出形参为函数的输出参数。当输出形参多于1个时，则应该用方括号括起来。,9,函数调用,函数调用的一般格式是：输出实参1,实参2,.=函数名(输入实参1,实参2,.）例：k,xk,yk,p=jhnewtonqx(1.5,20),10,文件保存位置,设置当前目录位置 cd C:Matlab7work%设置Matlab当前目录为：C:Matlab7work,11,Matlab基础（实验二）,1.矩阵输入2.特殊矩阵3.误差分析4.内置函数,12,矩阵输入,a=%建立空矩阵a,可在workspace中编辑 a=1,2,3;4 5 6%分行“,”或空格 分列“;”a=1:5%建立15矩阵(向量)a=(1,2,3,4,5)a=0:pi/5:2*pi%(pi=)建立111矩阵（向量）a：初始元0，步长/5,终止元2【注意“：”的用法】,13,常用矩阵,1、零矩阵 z=zeros(3,4)%产生 34零矩阵z；zeros(5)%产生 5阶零矩阵；2、单位矩阵 E=eye(6)%产生6阶单位矩阵E；2、幺矩阵 ones(3,2);%产生3行2列元素为1的矩阵；3、随机矩阵 r=rand(7)%产生01间分布的随机矩阵r s=round(rand(7)*30)%产生7阶030间均匀分部的随机矩阵s,14,特殊矩阵,1、Hilbert矩阵 h=hilb(5)%产生 5阶Hilbert矩阵h h=sym(hilb(5)%产生 5阶Hilbert符号矩阵h2、Vandermonde矩阵 v=vander(1,2,3,4,5)%产生5阶Vandermonde矩阵v3、魔方矩阵 m=magic(3)%产生3阶魔方矩阵3、Toeplitz矩阵 t=toeplitz(0:-1:-5,0:5),15,矩阵运算,a=magic(3);b=round(10*rand(3);r=(1:3);a+b%矩阵加法 a-b%矩阵减法 a*b%矩阵乘法 n=inv(a)%inv(a)=a(-1)=a的逆 a/b%a/b=a*inv(b)=a*b(-1)x=ar,a*x%x=ar=inv(a)*r=a(-1)*r c=a%a=a的（共轭）转置,16,元素运算,a=magic(3);b=round(10*rand(3);r=(1:3);a.*b%矩阵对应元素乘积 a b a.*b,a*b%比较 a.*b,a*b a./b%矩阵对应元素右除 a.b%矩阵对应元素左除 a.(1/3)%矩阵对应元素的立方根,17,矩阵索引,a=1:30%产生120行向量 a=reshape(a,5,6)%变更a的结构为35的矩阵 a(3,2)%取元素a（3,2）a(3,:)%取a第3行，取a第二列 a(:,2)%取a第2列 b=a(1,3,3:5)%取a的1,2行，2,4,5列元 m,n=size(a)%输出a的行列数m、n b=a(end,:)%b取a的最后一行,18,矩阵操作,a=1:36;a=reshape(a,6,6);t=toeplitz(0:-1:-5,0:5);c=diag(a)%提取a的列向量 b=diag(c)%b为以c为对角元的对角阵 d=diag(diag(a)d=diag(diag(t,-1),-1)l=tril(a)%l为a的下三角阵 u=triu(a,1)%u为a从第1条对角线的开始取的上三角阵 l=tril(a,-2)%l为a从第-2条对角线的开始取的下三角阵,19,矩阵分析,rank(a)%a的秩 norm(a)%a的2-范数 norm(a,inf)%a的无穷范数 cond(a)%a的谱（2-）条件数 cond(a,1)%a的1-条件数,20,矩阵分解与特征值,L,U=lu(a)%矩阵a的LU分解 Q,R=qr(a)%矩阵a的QR分解 eig(a)%矩阵a的特征值向量 V,D=eig(a)%D主对角线元素为A的全部近似特征值%V第k列元素为对应于A的特征值D(k,k)的特征向量,21,最大元搜索,m=magic(5)x,q=max(max(m)%矩阵m每列最大元素向量构成x，x向量最大元下标 C,p=max(m(:,q)%向量m(:,q)最大元素C，m(:,q)最大元下标p%例：搜索A的非对角线上绝对值最大元素A(p,q)A=abs(m-triu(m);x,q=max(max(A);.M,p=max(A(:,q);p,q,A(p,q),22,条件搜索,m=magic(5),n,n=size(m);s=find(m24,1)%搜索m中大于24的元素的前1个，返回此元素下标s q=ceil(s/n)%计算第s个元素所在列 p=s-n*(q-1)%计算第s个元素所在行,23,二维绘图,（一）plot 最基本的二维图形指令plot的功能：plot命令自动打开一个图形窗口Figure 用直线连接相邻两数据点来绘制图形根据图形坐标大小自动缩扩坐标轴，将数据标尺及单位标注自动加到两个坐标轴上，可自定坐标轴，可把x,y 轴用对数坐标表示,24,如果已经存在一个图形窗口，plot命令则清除当前图形，绘制新图形可单窗口单曲线绘图；可单窗口多曲线绘图；可单窗口多曲线分图绘图；可多窗口绘图可任意设定曲线颜色和线型可给图形加坐标网线和图形加注功能,25,plot的调用格式,plot(x)缺省自变量绘图格式，x为向量,以x元素值为纵坐标，以相应元素下标为横坐标绘图 plot(x,y)基本格式，以y(x)的函数关系作出直角坐标图，如果y为nm的矩阵，则以x 为自变量，作出m条曲线plot(x1,y1,x2,y2)多条曲线绘图格式,26,plot(x,y,s)开关格式，开关量字符串s设定曲线颜色和绘图方式，使用颜色字符串的前13个字母，如 yellowyel表示等。或plot(x1,y1,s1,x2,y2,s2，),27,S的标准设定值,字母 颜色 标点 线型 y 黄色 点线 m 粉红 圈线 c 亮蓝 线 r 大红 字线 g 绿色 实线 b 蓝色 星形线 w 白色：虚线 k 黑色(-)点划线,28,+、o、*、.、x、v、s quare 正方形 d iamond 菱形 p entagram 五角星 h exagram 六角星,matlab线形,29,单窗口多曲线绘图,x=0:pi/100:2*pi;y=sin(x);y1=cos(x);plot(x,y,r:,x,y1,gp),30,子图绘制,x=-3:0.1:3;y1=x;y2=x.2;y3=x.3;y4=x.4;subplot(2,2,1);plot(x,y1);title(y1=x)subplot(2,2,2)plot(x,y2);title(y2=x2)subplot(2,2,3)plot(x,y3);title(y3=x3)subplot(2,2,4)plot(x,y4);title(y4=x4),31,图形标注,标题函数：tiltle x=0:0.1:5;y=exp(-0.2*x).*sin(x);plot(x,y)title(it e0.2xsin(x),FontWeight,Bold),32,图形标注,坐标轴标注：xlabel,ylabel x=1990:2:2000;y=1.25 0.81 2.16 2.73 0.06 0.55;xin=1990:0.2:2000;yin=spline(x,y,xin);plot(x,y,ob,xin,yin,-.r)title(1990年到2000年某地区年平均降水量图)xlabel(it 年份,FontSize,15)ylabel(降雨量,FontSize,8),33,图形标注,曲线标注：legend x=0:0.02*pi:2*pi;y1=sin(x);y2=cos(x);y3=sin(3*x).*cos(x);plot(x,y1;y2;y3)axis(0 2*pi-1.5 1.5)legend(sin(x),cos(x),sin(3x)cos(x),-1),34,TeX字符,bf%粗体it%斜体rm%正常字体 fontname%指定字体fontsize%指定字号x%上标_x%下标pi;lambda alpha;beta;,35,文本标注 text,x=linspace(-1,1,500);f=exp(x);l2=0.9963+1.1036*x+0.5367*x.*x;l3=0.9963+0.9979*x+0.5367*x.*x+0.1761*x.*x.*x;plot(x,f,-b,x,l3,-r);xlabel(0 leq x leq1);ylabel(y);title(bf ex的3次最佳平方逼近多项式,fontname,隶书,fontsize,12)legend(it f(x)=ex,it l_2=0.9963+0.9979*x+0.5367*x2+0.1761*x3,2)text(0.96,exp(0.96),it f(x)=ex,rightarrow,.HorizontalAlignment,right,fontsize,12,color,0,0,1,fontname,times new roman);text(0.6,0.9963+0.9979*0.6+0.5367*0.36+0.1761*0.36*0.6,it l_3=0.9963+0.9979*x+0.5367*x2+0.1761*x3,rightarrow,.HorizontalAlignment,right,fontsize,12,color,1,0,0,fontname,times new roman);,36,多项式函数,p=1:5;v=1:3;poly2sym(p)r=roots(p)poly(v)y=polyval(p,1)y=polyval(p,p)conv(p,v)q,r=deconv(p,v),37,多项式插值,p=polyfit(x,y,n)%将数据点x，y拟合为n次多项式p x=linspace(-5,5,11);y=5./(1+x.*x);p=polyfit(x,y,10);%10次插值多项式t=linspace(-5,5);z=polyval(p,t);plot(t,z)hold on;x=linspace(-5,5,5);y=5./(1+x.*x);p=polyfit(x,y,4);%4次插值多项式t=linspace(-5,5);z=polyval(p,t);plot(t,z,r),38,多项式插值,p=polyfit(x,y,n)%将数据点x，y拟合为n次多项式p x=linspace(-5,5,11);y=5./(1+x.*x);p=polyfit(x,y,10);t=linspace(-5,5);z=polyval(p,t);plot(t,z)yi=interp1(x,y,xi,method)%计算数据点(x,y)按method指定的插值函数在xi点的值，method参数：linear分段线性插值spline默认三次样条插值;cubic三次Hermite插值；x=0:10;y=sin(x);xi=0:0.25:10;yi=interp1(x,y,xi);%默认为一维线性插值plot(x,y,o,xi,yi),39,多项式插值,x=0:10;y=sin(x);xi=0:.25:10;yi=sin(xi);grid ony2=interp1(x,y,xi,spline);plot(x,y,xi,yi,xi,y2,-.)%三次样条插值y3=interp1(x,y,xi,cubic);%三次多项式插值plot(x,y,o,xi,yi,xi,y3,-.),40,多项式拟合函数,%计算给定数据点的n次拟合多项式：x=0.5,1,1.2,1.4,1.6,1.8,2 5;y=1.44./(x.*x)+0.24*x;dxsnh(x,y,3),41,最小二乘逼近多项式,%计算sin(2x)在区间 0,1 三次最佳平方逼近多项式，作出最佳平方逼近多项式图形。p,err,errinf=zjpfbj(3),42,匿名函数,f=(x)(x*exp(x)%创建以 x 为输入参数的匿名函数：xex y=feval(f,2)%计算函数：f(x)=xex 在 x=2点的函数值 f(2)=2e2,43,积分,F=int(x*exp(x)%求以 x 为积分变量的函数 xex 的原函数F I=int(x*exp(x),1,2)%求以 x 为积分变量的函数 xex 的在积分区间1，2上的定积分I,44,微分,syms x f=exp(1/x)%创建符号函数 f(x)=e(1/x)d=diff(f)subs(d,x,2),