4 数值计算实验 4.1 基础训练 1. 方程求根 编程调用 fzero 求解方程 在 附近的一个根，并将所求根赋给 变量 xp，并求出该方程所有的数值根赋给变量 r。 解： function [xp,r]=myfun fun=@(x)2*x.^3-3*x.^2+4*x-5; x = -5:0.1:10; plot(x,fun(x)) %绘制函数曲线 [xp,val,flag]=fzero(fun,2); r=roots([2,-3,4,-5]); 2. 求解二阶微分方程 编写函数调用 ode 工具箱函数返回 x 在点 0:0.1:5 处的函数值，用列向量 r 存储这些函 数值, 此列向量为 double 型数组. 解： function r=myfun x=0:0.1:5; [T,Y]=ode23(@dfun,0:0.1:5,[2;0]); plot(T,Y(:,1),'-') r=Y(:,1); function dy=dfun(t,y) %x--y(1), dx/dt--y(2) dy=[y(2);20*(1-y(1)^2)*y(2)+0.5*y(1)]; 3. 二次多项式拟合 某种产品在生产过程中的性能指标 y 与它所含的某种材料的含量 x 有关，现将试验所得 16 组数据记录列于下表。 X 22.09 20.05 26.5 36.23 10.46 38.2 Y X Y 24.13 2.75 40.27 26.24 3.53 42.27 28.11 11.67 44.07 30.29 29.98 46.05 445.1 32.09 52.26 48.47 34.23 87.19 50.08 552.84 631 128.11 176.24 235.17 300.25 365.66 要求拟合 y 与 x 的函数关系。用多项式拟合函数 polyfit 进行二次多项式拟合。编写函数 0543223xxx2x0)0(;2)0(05.0)1(20222xxxdtdxxdtxd 

文件返回 2 个参数： 第 1 个返回参数为二次多项式系数组成的行向量 p（元素由高次到低次排列）； 第 2 个返回参数为拟合函数在点 x=25:0.4:60 处的函数值（用 1 个行向量表示）。 程序文件第 1 行参考格式如下： function [p,v]= myfun 解： function [p,v]= myfun % 曲线拟合 % y= 2.00 + 1.00*(x-25.00).^2 % 多项式系数: 1 -50 627 % 多项式系数: % y=1.00 x^2 + -50.00 x + 627.00 x=[ 20.05 22.09 24.13 26.24 28.11 30.29 32.09 34.23 ... 36.23 38.2 40.27 42.27 44.07 46.05 48.47 50.08]; y=[ 26.5 10.46 2.75 3.53 11.67 29.98 52.26 87.19 ... 128.11 176.24 235.17 300.25 365.66 445.1 552.84 631]; plot(x,y,'o') p=polyfit(x,y,2) v=polyval(p, 25:0.4:60) %得到多项式在数据点处的值 4. 编写一个 function 文件求定积分 1) 求定积分 . 2) 求二重积分 要求该function文件返回2个参数[v1,v2]分别表示所得积分结果. 解： function [v1,v2]=testmain v1 = quad(@myfun,0,3) v2=dblquad(@fun,-1,1,-2,2) %定义函数 function y = myfun(x) y = exp(x)./(x.^3+2*exp(x)+1); function r=fun(x,y) r=(exp(y).*x.^2+exp(x).*y.^2).*(x.^2+y.^2/4<=1); 30312dxexexxdxdyyexeyxxy)(142222

4.2 实验任务 一．实验任务 请用 Euler 法和 Matlab 函数 ode23 求解下列微分方程： . 并将 Euler 求解结果与 Matlab 的 ode23 函数求解结果対比。 二. 实验目的 认识 Euler 法。熟悉 Matlab 解微分方程数值解的函数. 三. 实验过程 用差商近似微商，可得 从而将原问题化为如下形式： , . 可以采用不断迭代的方式估算不同时刻 y 的近似值。 参考程序如下： function testmain % Euler法 解常微分方程 clear x0 = 0; y0 = 10; h = 0.01; x= 0:h:300; y= zeros(size(x)); % x(i) <--> y(i) y(1) = y0; for i=1:length(x)-1 % Euler方法解微分方程 y(i+1) = y(i)+h*0.02*(1-0.001*y(i))*y(i); end plot(x(1:5:end),y(1:5:end),'--') [T,Y]=ode23(@dfun,[0 300],10); %调用ODE函数 hold on plot(T(1:2:end),Y(1:2:end),'r*') 10)0()001.01(02.0yyydtdy10)0()001.01(02.0yyydtdyttyttydtdy)()(ttytytytty)())(001.01(02.0)()(

function dy=dfun(t,y) dy=0.02*(1-0.001*y).*y; 四. 实验自评与改进方向 1. 总体完成较为顺利。 2. 可以绘制动态变化过程，并制作动画视频。 五. 实验体会，收获及建议 1. 可以举一个实际问题中的微分方程或微分方程组的例子。