常州大学实习报告
2 详细求解
2.1必做题
2.1.2 求解常微分方程x’’=-x’+x+1，x’(0)=1，x(0)=0（使用函数dsolve）
>>dsolve('D2x=-Dx+x+1','Dx(0)=1','x(0)=0') ans = exp(t*(5^(1/2)/2 - 1/2))*((3*5^(1/2))/10 (5^(1/2)*exp(-t*(5^(1/2)/2 + 1/2))*(5^(1/2) - 3))/10 - 1
+ 1/2) +
2.1.3 已知t=an2+bn，测得对应数据如下：（多项式插值interp1） t=[0,20,40,60,80,100,120,140,160,183.5];
n=[0,1153,2045,2800,3466,4068,4621,5135,5619,6152]; 试求a和b的值。 >>t0=[0,20,40,60,80,100,120,140,160,183.5]; n0=[0,1153,2045,2800,3466,4068,4621,5135,5619,6152]; n=0:0.001:6152; t=interp1(n0,t0,n,'spline'); p=polyfit(n,t,2) p = 0.0000 0.0144 0.0631 2.1.4 请用梯形法、辛普森法分别计算积分值 （trapz、quad） 2.1.5 计算二重积分?1200
?(x2?y2?xy?2x?y?1)dxdy （使用函数dblquad）
>>dblquad('x.^2+y.^2+x*y+2*x+y+1', 0, 1, 0, 2) ans = 10.3333 >>f=inline('sqrt(x.^2+x+1)','x'); >> quad(f,0,1) ans = 1.3369 >> x=0:0.01:1;y=sqrt(x.^2+x+1);trapz(x,y) ans = 1.3369 第 2 页 共 16 页
常州大学实习报告
2.1.6 矩阵M=[1,2,6; 4,2,7; 8,9,3]，求M的LU分解，QR分解，对角阵，特征值分解。（使用函数lu、qr、shol、eig）
2.1.7 a=3,A=4,b=a2,B=b2-1,c=a+A-2B,C=a+B+2c，求C。（使用函数solve） >>M=[1,2,6; 4,2,7; 8,9,3]; lu(M) ans = 8.0000 9.0000 3.0000 0.5000 -2.5000 5.5000 0.1250 -0.3500 7.5500 >>qr(M) ans = 9.0000 9.1111 6.4444 -0.5000 2.4470 -2.3360 -1.0000 -0.4719 6.8566 >>[D,X]=eig(M) D = -0.4111 -0.7719 -0.4992 -0.5484 0.6237 -0.3914 -0.7282 0.1229 0.7730 X = 14.2944 0 0 0 -1.5712 0 0 0 -6.7232 >> syms a b c A B C eq1='a=3'; eq2='A=4'; eq3='b=a^2'; eq4='B=b^2-1'; eq5='c=a+A-2*B'; eq6='C=a+B+2*c'; x=solve(eq1,eq2,eq3,eq4,eq5,eq6); >> x.C ans = -223 >>diag(M) ans = 1 2 3 2.1.8用两种方法求解Ax=b的解。（A为四阶随机矩阵，b为四阶向量，自己创建）。 方法一： >> A=rand(4) A = 0.9501 0.8913 0.8214 0.9218 0.2311 0.7621 0.4447 0.7382 0.6068 0.4565 0.6154 0.1763 0.4860 0.0185 0.7919 0.4057
第 3 页 共 16 页
>> B=[1;2;3;4];x=inv(A)*B x = -6.4970 1.9163 10.9584 -3.8362 常州大学实习报告
方法二：
>> A(:,5)=B;rref(A) ans = 1.0000 0 0 0 -6.4970 0 1.0000 0 0 1.9163 0 0 1.0000 0 10.9584 0 0 0 1.0000 -3.8362
12.1.9 f(x)?x3?(x?0.8)2/(x?1.25)3?5(x?)，用两种方法求函数的根，并求
x其极值与零点。 方法一：
>> solve('x^3+(x-0.8)^2/(x+1.25)^3-5*x-5/x');x=double(ans); >> for i=1:length(x) xx(i)=isreal(x(i));
end
>> x=x(xx) x =
2.4156
方法二：
>> fzero('x.^3+(x-0.8).^2/(x+1.25).^3-5*(x+1/x)',3) ans = 2.4156
求极值：
>> syms x y >> y=x.^3+(x-0.8).^2/(x+1.25).^3-5*(x+1/x); >> dydx=diff(y); >> solve(dydx);x=double(ans); >> for i=1:length(x) xx(i)=isreal(x(i)); end >> t=x(xx) t = -0.4694 -2.4039
>> x=t(1);y1=subs(y) y1 = 16.2832 >> x=t(2);y2=subs(y) y2 = -6.4732 >> z=diff(diff(y)); >> x=t(1);z1=subs(z) z1 = 205.8164 >> x=t(2);z2=subs(z) z2 = -53.5382 第 4 页 共 16 页
常州大学实习报告
所以函数f(x)存在一个极大值点x1=-2.4039,极大值为-6.4732；一个极小值点x2=-2.4039，极小值为16.2832
?x2?x?6,x?0且x??4?22.1.10 f(x)的定义如下：f(x)??x?5x?6,0?x?10,x?2且x?3，写一个matlab
2?x?x?1,其它?函数func1实现该函数。
2.1.11 写一个MATLAB小程序，求出最小的 n 值，使得 n! > realmax。
function func1 (x) if x<0&&x==-4 y=x^2+x-6; elseif x>=0&x<10&x~=2&x~=3 y=x^2-5*x+6; else y=x^2-x-1; end fprintf('%d\\n',y); function f=realmax() s=1;n=1; while(1) n=n+1; s=s*n; if s>realmax break; end end fprintf('最小的整数为：%d\\n',n)
2.1.12 写一个 MATLAB 函式 myfun.m 来计算下列方程式：
y = 0.5*exp(x/3)-x*x*sin(x)
其中 x 是函式的输入，y 是函式的输出。你的函式必须能够处理当 x 是纯量或是向量的两种情况。
function f=myfun(x) n=length(x); for i=1:n f(i)=0.5*exp(x(i)/3)-x(i)*x(i)*sin(x(i)); end
2.1.13 写一个 MATLAB 函式 pifun.m 来计算下列级数：
f(n) = 4*(1 - 1/3 + 1/5 - 1/7 + 1/9 - 1/11 + ...)
其中 n 为函式的输入，级数和 f(n) 则是函式的输出。此外，你必须进行下列事项：
（1） 使用 tic 和 toc 指令来测量 pifun(100000) 的计算时间。如果你不知道如何使用这两个指令，请使用 help tic 及 help toc 来查出它们的用法。
第 5 页 共 16 页
常州大学实习报告
我的计算机是 Pentium-450，所得的计算时间约为 2 秒。请说明你的计算机规格以及其计算时间。
（2）使用 flops 指令来测量 pifun(100000) 所用到 floating point operations 的计算次数。如果你不知道如何使用这个指令，请使用 help flops 来查出它的用法。
function [sum,time,flop]=pifun(n) %f(n) = 4*(1 - 1/3 + 1/5 - 1/7 + 1/9 - 1/11 + ...) tic sum=0; i=0; flag=1; for i=0:1:n sum = sum + 4*(flag/(2*i+1)); flag = -flag; end time=toc 2.1.15 求下列函数的极小点：
222 1)f?X??x1?4x2?9x3?2x1?18x2
>> f='x(1)^2+4*x(2)^2+9*x(3)^2-2*x(1)+18*x(2) '; [x,fm]=fminsearch(f,[0,0,0]) x = 1.0000 -2.2500 0.0000 fm = -21.2500 3222) f?X??x1?x2?2x1x2?x1?2x2
2
>> f='x(1)^2+3/2*x(2)^2-2*x(1)*x(2)+x(1)-2*x(2) '; [x,fm]=fminsearch(f,[0,0]) x = 0.499977193672681 1.000008683829168 fm = -0.749999998970666
第 6 页 共 16 页