科学计算与编程.docx
- 文档编号:11184399
- 上传时间:2023-02-25
- 格式:DOCX
- 页数:17
- 大小:18.56KB
科学计算与编程.docx
《科学计算与编程.docx》由会员分享,可在线阅读,更多相关《科学计算与编程.docx(17页珍藏版)》请在冰豆网上搜索。
科学计算与编程
1.求解非线性方程的所有解:
fzero;
题:
求解y=x.*sin(1./x)-0.2*exp(-x)在[0,2]的解。
1.写出函数:
functiony=myfun01(x)
y=x.*sin(1./x)-0.2*exp(-x);
end
2.再新建一个文件:
clc
closeall
clearall
figure
fplot(@myfun01,[0,2])
gridon(画图)
x1=fzero(@myfun01,[0.2,0.4])(找出零点所在的区间)
myfun01([x1])(验算是否是方程的解)
2.使用弦截法求解非线性方程的解;
题3-1为例
步骤:
1.用function写出求解函数;
functiony=myfunction(x)
y=0.5*exp(x/3)-sin(x);
end
2.用function创建弦截法求解非线性方程的程序secant函数:
secant函数如下:
functionr=secant(f,a,b,n)
fa=f(a);
fb=f(b);
fork=1:
n
m=a-(b-a)*fa/(fb-fa);
fm=f(m);
if(fm==0)
break
elseif(fa*fm<0)
b=m;fb=fm;
else
a=m;fa=fm;
end
end
r=a-(b-a)*fa/(fb-fa);
3.在主文件下(新建一个文件做主文件)画图确定方程解的个数及解所在区间
clc
closeall
clearall
f=@myfunction;
figure
(1)
fplot(f,[0,10])
gridon
x1=secant(f,0,1,100)
x2=secant(f,1,2,100)
3.使用fsolve求解非线性方程组的解;
P75习题3-4
程序如下:
1、先建一个myfun文件
functionF=myfun01(X)
x=X(01);
y=X(02);
F
(1)=x.*exp(x*y+0.8)+exp(y^2)-3;
F
(2)=x^2-y^2-0.5*exp(x*y);
2、在建Edit文件编写主程序
clc
clearall
closeall
[fsolve(@myfun01,[11])
fsolve(@myfun01,[1-1])
fsolve(@myfun01,[-11])]求出三个区间的三个解
4.求解线性方程组的解(恰定方程组、超定方程组);
第四题:
(4-1到4-6都可以用这个方法写出来)
本题为4-4:
clc
clearall
closeall
A=[1,1,2;-1,2,1;1,-1,2];
b=[4,4,-1]';
x=pinv(A)*b
5.使用多项式进行曲线拟合;
P101例1
拟合一次多项式:
clc
clearall
closeall
x=[1,1.2,1.4,1.6,1.8];
y=[1,1.0954,1.1832,1.2649,1.3416];
plot(x,y,’r.’,’markersize’,25)
gridon先画图,再接着写下面的
x=x’;
A=[X.^1;X.^0]’;
b=y’;
a=pin(A)*b或a=A\b
holdon写到这运行程序,第一个数是x的系数,第二个是常数项
fplot(‘0.4264*x+0.5801’,[0.8,2])图线区间
拟合二次多项式:
第5.4题
1、程序:
clc
closeall
clearall
x=[0.51.01.52.02.53.0];
y=[1.752.453.814.808.008.60];
figure
plot(x,y,'r.','markersize',25)
gridon
A=[x.^2;x.^1;x.^0]';
b=y';
a=A\b
holdon先写到这里,运行,得到三个数,第一个是二次项系数,以此类推,写出下式。
fplot('0.4900*x^2+1.2501*x+0.8560',[0.53])
运行得到连线图。
6.使用经验公式进行曲线拟合:
a幂函数:
y=βxα
以题5-5为例
程序:
clc
closeall
clearall
x=[2.02.22.42.62.83.03.23.43.63.84.0];
y=[40.154.373.198.8133.4180.2243.0328.0442.9597.7806.8];
figure
(2)
plot(x,y,'.r','markersize',30)
gridon
holdon
X=log(x);求y=βxα中β和α的值
Y=log(y);
a=polyfit(X,Y,1)
a1=a
(1);a0=a
(2);
alpha=a1α的值
beta=exp(a0)β的值
f='1.6950*x.^(4.3381)';(根据结果写出函数)
fplot(f,[2,4.5])(画出函数图像)
b指数函数:
y=αeβx
以题5-2为例
程序:
clc
closeall
clearall
x=[0.01290.02470.05300.15500.30100.47100.80201.27001.43002.4600];
y=[9.56008.18455.26162.79172.26111.73401.23701.06741.11710.7620];
figure
(1)
plot(x,y,'.r','markersize',30)
gridon
holdon
X=x;求y=αeβx中β和α的值
Y=log(y);
a=polyfit(X,Y,1)
a1=a
(1);a0=a
(2);
alpha=exp(a0)α的值
beta=a1β的值
f='4.4717*exp(-0.9238*x)';(根据结果写出函数)
fplot(f,[0,3])(画出函数图像)
c其他函数:
y=a+bx2等
以P1346-1为例
求y=a+bx
clc
closeall
clearall
x=[0.10.40.50.70.80.9];
y=[0.610.920.991.521.472.03];
A=[x.^1;x.^0]';(若求y=a+bx2则改为A=[x.^2;x.^0]'若求y=a0+a1x+a2x2则改为A=[x.^2;x.^1;x.^0]')
plot(x,y,'r.','markersize',30)
gridon
b=y';
a=pinv(A)*b
holdon
fplot('1.6577*x+0.3173',[0.1,1])
以P134【5-4】为例求y=a0+a1x+a2x2
clc
closeall
clearall
x=[0.51.01.52.02.53.0];
y=[1.752.453.814.808.008.60];
A=[x.^2;x.^1;x.^0]';
plot(x,y,'r.','markersize',30)
gridon
b=y';
a=pinv(A)*b
holdon
fplot('0.49*x^2+1.2501*x+0.8560',[0.5,3])
7.分段线性插值:
interp1(x,y,x0)
P160以6-2为例:
程序如下:
clc
clearall
closeall
x=[1,1.4,1.8,2.3,2.6,3.0];
y=[1.000,1.1832,1.3416,1.5166,1.6125,1.7321];
figure
(2)
plot(x,y,'.r','markersize',35)
gridon
holdon
interp1(x,y,1.69)(求分段线性插值近似f(1.69)=?
)
x0=1:
0.1:
3;(绘制插值函数的图像)
y0=interp1(x,y,x0);
plot(x0,y0,'g')
8.多项式插值(拉格朗日插值):
log_polyfit(x,y,x0),polyfit/polyval
P160以6-1为例:
1.先建一个拉格朗日函数:
functiony_fit=lag_polyfit(x,y,x_fit)
if(length(x)==length(y))
n=length(x);
else
disp('error');
return;
end
y_fit=0;
fori=1:
n
item=1;
forj=1:
i-1
item=item.*(x_fit-x(j))/(x(i)-x(j));
end;
forj=i+1:
n
item=item.*(x_fit-x(j))/(x(i)-x(j));
end;
si=item.*y(i);
y_fit=y_fit+si;
end
2.再建主函数
clc
clearall
closeall
x=[1.0,1.6,2.2,2.8,3.4,4];
y=[0.8415,0.9996,0.8085,0.3350,-0.2555,-0.7568];
figure
plot(x,y,'.r','markersize',35)
gridon
holdon
lag_polyfit(x,y,1.8)(求拉格朗日插值近似f(1.8)=?
)
x0=1:
0.1:
4;
y0=lag_polyfit(x,y,x0);
plot(x0,y0,'g')
a=polyfit(x,y,length(x)-1)(求拟合多项式的系数)
f='-0.0062*x^5+0.0998*x^4-0.4780*x^3+0.5117*x^2+0.5736*x+0.1406'(根据a的结果写出函数,第一个数是最高次的系数依此类推)
fplot(f,[1,4])(画出插值函数图)
9.三次样条插值:
spline(x,y,x0)
P160习题6-3
程序如下:
clc
clearall
closeall
x=[1.1,1.3,1.7,2.2,2.5];
y=[1.2100,2.1832,2.3416,2.4156,2.7155];
y_fit=spline(x,y,1.6);
disp(y_fit)(用三次样条插值法求点的值f(1.6))
figure
(1)
plot(x,y,'.r','markersize',25);
x_fit=1.1:
0.01:
2.5;
y_fit=spline(x,y,x_fit);
plot(x_fit,y_fit,'linewidth',4);
gridon(绘制插值函数的图像)
得到的答案:
f(1.6)=2.3935
10.解析法求积分:
int
a不定积分
b定积分
c多重定积分
P1857-1题
(1)不定积分
clc
clearall
closeall
D1=int('(1/x^2)*cos(1/x)','x')
没有vpa
(2)定积分
clc
clearall
closeall
D2=int('x/(1+cos(2*x))','x',0,pi/4)
vpa(D2,8)
(3)多重定积分(3重定积分)
clc
clearall
closeall
D3=int('1/(1+exp(x))^2','x',0,inf)
vpa(D3,8)
(4)4重定积分
clc
clearall
closeall
D4=int('1/(x^2+4*x+9)','x',-inf,inf)
vpa(D4,8)
7-6题二重定积分:
clc
clearall
closeall
D5=int('exp(-x^2/2)*sin(x^2+y)','x',-2,2)
D6=int(D5,'y',-1,1)
vpa(D6,8)
7-7题
clc
clearall
closeall
D7=int('4*x*z*exp(-x^2*y-z^2)','x',0,2)
D8=int(D7,'y',0,pi)
D9=int(D8,'z',0,pi)
vpa(D9,8)
11.数值法求定积分:
a梯形法、辛普森法、复合梯形法
例题7-2(梯形法):
clc
clearall
closeall
a=0;
b=pi/4;
fa=a/(1+cos(2*a));
fb=b/(1+cos(2*b));
h=b-a;
I=h/2*(fa+fb)
II=int('x/(1+cos(2*x))','x',0,pi/4)
vpa(II)
例题7-3(复合梯形法):
clc
clearall
closeall
a=0;
b=pi/4;
n=500;
h=(b-a)/n;
x=a:
h:
b;
f=x./(1+cos(2*x));(此处需点乘)
I=h/2*(f
(1)+f(n+1))
II=h*sum(f)-I
例题7-4(辛普森法):
clc
clearall
closeall
formatlong
a=0;
b=pi/4;
x0=a;
x1=(a+b)/2;
x2=b;
f0=x0/(1+cos(2*x0));
f1=x1/(1+cos(2*x1));
f2=x2/(1+cos(2*x2));
I=1/6*(b-a)*(f0+4*f1+f2)
II=int('x/(1+cos(2*x))','x',0,pi/4)
b三次样条法:
spline/fnint/fnval
以【7-8】为例
clc
clearall
closeall
formatlong
x=0:
0.6:
3.0;
y=[0,0.5646,0.9320,0.9738,0.6755,0.1411];
s=spline(x,y);
F=fnint(s);
I=fnval(F,3)-fnval(F,0)
c多重定积分数值函数
P185以7-6为例
D1=int('exp(-x^2/2)*sin(x^2+y)','x',-2,2)
D2=int(D1,'y',-1,1)
vpa(D2)(解析法)
dblquad('exp(-x.^2/2).*sin(x.^2+y)',-2,2,-1,1)(数值法)
以【7-7】为例
D1=int('4*x*z*exp(-x^2*y-z^2)','x',0,2);
D2=int(D1,'y',0,pi);
D3=int(D2,'z',0,pi);
vpa(D3)(解析法)
triplequad('4*x.*z.*exp(-x.^2.*y-z.^2)',0,2,0,pi,0,pi)(数值法)
12、解析法求微分:
a)1阶导函数、2阶导函数、3阶导函数:
diff【8-1】
b)1阶导数值、2阶导数值、3阶导数值:
diff/subs【8-1】
程序:
clc
closeall
clearall
df1=diff('exp(x)*cos(x)','x',1)(f(x)对x的一阶导数)
df2=diff('exp(x)*cos(x)','x',2)(f(x)对x的二阶导数)
df3=diff('exp(x)*cos(x)','x',3)(f(x)对x的三阶导数)
d1=subs(df1,'x','0.5')(当x=0.5时f(x)的一阶导数值)
d2=subs(df2,'x','0.5')(当x=0.5时f(x)的二阶导数值)
d3=subs(df3,'x','0.5')(当x=0.5时f(x)的三阶导数值)
df4=diff('exp(x)*cos(y)','x',2)(z对x的二阶导数)
df5=diff(df4,'y',1)(再对y求一阶导数)
13、数值法求微分:
A)三次样条法:
spline/fnder/fnval【8-4】
B)拉格朗日法:
polyfit/ployder/polyval【8-5】
a)三次样条法
x=3:
0.5:
6;
y=[1.7321,1.8708,2.0000,2.1213,2.2361,2.3452,2.4495];
s=spline(x,y);
dsn1=fnder(s,1);一阶导数
dsn2=fnder(s,2);二阶导数,若要求三阶就再加上dsn3=fnder(s,3),同理继续
fnval(dsn1,4)x=4处
fnval(dsn2,4)x=4处,三阶将dsn2改为dsn3
b)拉格朗日法
clc
x=3:
0.5:
6;
y=[1.7321,1.8708,2.0000,2.1213,2.2361,2.3452,2.4495];
L=polyfit(x,y,length(x)-1);
dsn1=polyder(L);
dsn2=polyder(dsn1);
polyval(dsn1,4)在x=4处的一阶导数
polyval(dsn2,4)
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 科学 计算 编程