h080504.m

(有源代码)数值分析作业,本文主要包括两个部分,第一部分是常微分方程(ODE)的三个实验题,第二部分是有关的拓展讨论,包括高阶常微分的求解和边值问题的求解(BVP).文中的算法和算例都是基于Matla

%-- Unknown date --%
y=x;
o=ones(1,nkx);
kx=x' * o;
ky=(y' *o)';
kr=abs(complex(kx,ky));
kf= ((rhos * h + i * rho ./sqrt(k^2-complex(kr,0.0).^2))* omega^2/D).^0.25;
kfa=kf*a;
ss=-rho*omega^2./(D*(2*pi)^2 *(complex(kr,0).^4 -kfa.^4));
wavei(dba(ss)',x,y)
shading('flat')
wavei(dba(ss)',x,y)
surf (ss); figure(gcf)
surf(dba(ss))
shading('flat')
colormap
colormap('jet')
colormap jet
colormap cool
colormap jet
wavei(dba(ss)',x,y)
open dba
open wavei
isreal(ss)
isreal(dba(ss))
isreal(dba(ss)')
which
whos
cc
whats
what
circ
sharm
cc
monopole
edit monopole
dipole
hard_sphere
soft_sphere
helical
open helical
cylrad
open cylrad
cc
syms u
syms B D
f0=(sinh(B*u)/B -sin(D *u) /D )/(B^2 + D^2)
mhelp sinh
clc
help prettify
prettify(f0)
pretty(f0)
fi
cd
fi
diff(f0,u)
diff(f0,u,2)
fi(2)
y=fi(2)
fi(2)
fi(0)
fi(1)
fi(2)
subs(ans,u,0.3)
vpa(ans)
doc sub
cc
doc subs
cc
%-- 08-4-2  下午10:15 --%
help pdepe
%-- 08-4-3 上午9:45 --%
doc pdepe
pdex4
whos
pack
cc
pdex1
plot (x, 'DisplayName', 'x', 'YDataSource', 'x'); figure(gcf)
plot (z, 'DisplayName', 'z', 'YDataSource', 'z'); figure(gcf)
cc
plot (F, 'DisplayName', 'F', 'YDataSource', 'F'); figure(gcf)
plot (DfDy, 'DisplayName','DfDy', 'YDataSource', 'DfDy'); figure(gcf)
cc
plot (M0, 'DisplayName','M0', 'YDataSource', 'M0'); figure(gcf)
plot (UM, 'DisplayName','UM', 'YDataSource', 'UM'); figure(gcf)
plot (VM, 'DisplayName','VM', 'YDataSource', 'VM'); figure(gcf)
plot (J, 'DisplayName','J', 'YDataSource', 'J'); figure(gcf)
plot (AlgVar, 'DisplayName', 'AlgVar', 'YDataSource', 'AlgVar'); figure(gcf)
DFDY(AlgVar,AlgVar)
J = DFDY(AlgVar,AlgVar)
plot (J, 'DisplayName','J', 'YDataSource', 'J'); figure(gcf)
eps*nnz(J)*condest(J)
nnz(J)
condest(J)
AlgVar
plot (x, 'DisplayName', 'x', 'YDataSource', 'x'); figure(gcf)
plot (z, 'DisplayName', 'z', 'YDataSource', 'z'); figure(gcf)
%-- 08-4-3 下午1:07 --%
%-- 08-4-5 下午9:20 --%
doc lsqcurvefit
cc
whos
size(myfun(x,xdata))
whos
plot(ydata1)
plot(xdata,ydata1)
plot (xdata, 'DisplayName', 'xdata', 'YDataSource', 'xdata'); figure(gcf)
plot (ydata, 'DisplayName', 'ydata', 'YDataSource', 'ydata'); figure(gcf)
k = 1:10;
2*k
%-- 08-4-7 下午8:13 --%
help powergui
help dde
help dee
dee
%-- 08-4-8 下午7:03 --%
open untitled.mdl
cc
100*10^(-6)
format lont
format long
ans
format
ans
0.0001
cc
%-- 08-4-16 下午4:23 --%
look for van
lookfor van
doc vdpode
%-- 08-4-18 下午1:19 --%
ode45('eg63fun',[0 30],[1 0.5]);
[t,x]=ode45('eg63fun',[0 30],[1 0.5]);
plot(x(:,1),x(:,2))
open fuzzypid
cc
0.03sin(2*pi*3*t)
cc
Afp
cc
fuzzy
%-- 08-4-25 下午9:54 --%
s = signal('colpitts.dat','ascii')
view(s)
a = amutual(s,32);
help amutual
disp amutual
open amutual
a = amutual(s,32);
clear
clc
clear
clc
%-- 08-4-25  下午10:04 --%
s = signal('colpitts.dat','ascii')
view(s)
a = amutual(s,32);
view(a)
c = cao(s,8,4,3,1000);
%-- 08-4-25  下午10:11 --%
s = signal('colpitts.dat','ascii')
view(s)
a = amutual(s,32);
view(a)
c = cao(s,8,4,3,1000);
view(c)
e = embed(s,3,4)
view(e)
view(corrsum(e,-1, 0.05, 40 ,32));
takens_estimator(e, -1, 0.05 , 40)
view(largelyap(e,1000,300,40,2));
tstool
x = 1:100;
y = floor((x+rand(1,100))/pi);
xx = 1: 0.1:100;
yy = akimaspline(x,y,xx);
tsdemo
tsdemo2
tsinfo
cc
%-- 08-4-28  下午10:04 --%
open('函数处理的结果.fig')
%-- 08-4-29 下午9:51 --%
help laplace
solve('(sin(x)-x)/(1-cos(x)-0.5*x^2=tan(m)',m)
help solve
solve('(sin(a)-a)/(1-cos(a)-0.5*a^2=tan(x)')
solve('(sin(a)-a)/(1-cos(a)-0.5*a^2-tan(x)')
solve('(sin(a)-a)/(1-cos(a)-0.5*a^2-tan(x))')
solve('(sin(a)-a)/(1-cos(a)-0.5*a^2=tan(x))')
solve('(sin(a)-a)/(1-cos(a)-0.5*a^2)=tan(x)')
cc
%-- 08-4-30 下午3:00 --%
cd E:\matlabtemp
cc
a=3;
func = inline(' 1 + (y - x) .* a ');
func = inline(' 1 + (y - x) .* a ')
[x,y]=nm_ch_0501(-0.5,0.001);
plot(x,y)
size(x)
size(y)
[x,y]=nm_ch_0501(-0.5,0.001);
plot(x,y)
hold on
plot(x,exp(-0.5*x) +x)
plot(x,exp(-0.5*x) +x,'*r')
[x,y,r]=nm_ch_0501(-0.5,0.001);
plot(x,y,x,r)
cc
[x,y,r]=nm_ch_0501(1,0.1);
plot(x,y,x,r)
plot(x,y,x,r,'+r')
[x,y,r]=nm_ch_0501(10,0.1);
plot(x,y,x,r,'+r')
[x,y,r]=nm_ch_0501(50,0.1);
plot(x,y,x,r,'+r')
[x,y,r]=nm_ch_0501(50,0.05);
plot(x,y,x,r,'+r')
[x,y,r]=nm_ch_0501(50,0.005);
plot(x,y,x,r,'+r')
hold on
plot(x,abs(y-r))
plot(x,abs(y-r),'y<')
cc
cd E:\MATLAB7\Numerical Methods\Chapter5
cd 'E:\MATLAB7\Numerical Methods\Chapter5'
contents
help contents
type RK
cc
cd_temp
fix(3.5)
floor(3.4)
fix(-3.4)
floor(-3.4)
round(-3.4)
round(3.4)
maple('with[inttrans];')
maple('with[inttrans]:')
maple('with[inttrans]')
maple('laplace(diff(sin(3*t),t$2),t,s)')
maple('laplace(diff(f(3*t),t$2),t,s)')
maple('diff(f(t),t$2)')
maple('laplace(diff(f(t),`$`(t,2))),t,s)')
maple('laplace(diff(f(t),`$`(t,2)),t,s)')
maple('subs(laplace(f(t),t,s)=F(s),laplace(diff(f(t),`$`(t,2)),t,s))')
syms s
syms t
f(t)
syms f
f(t)
syms f(t)
sym('f(t)')
diff(f(t),2)
diff(ans,2)
f1=sym('f(t)')
diff(f1,2)
laplace(ans)
s1=laplace(f1)
s1=laplace(diff(f1,2))
subs(s1,laplace(f(t),t,s),X(s))
syms('X(s)')
X=sym('X(s)')
subs(s1,laplace(f(t),t,s),X)
subs(s1,laplace(f(t),t,s),1)
subs(f1,f,1)
syms s
subs(s1,laplace(f(t),t,s),1)
syms t s
subs(s1,laplace(f(t),t,s),1)
cc
feval(funtest1,1,2,0.5)
feval('funtest1',1,2,0.5)
[x,y]=rk_wjl('funtest1',0,1,1,0.01,-0.5)
plot(x,y)
[x,y]=rk_wjl('funtest1',0,1,1,0.01,-0.5)
[x,y]=rk_wjl('funtest1',0,1,1,0.01,-0.5);
plot(x,y)
hold on
plot(x,exp(-0.5*x) +x,'*r')
[x,y]=adams_wjl('funtest1',0,1,1,0.01,-0.5);
b
[x,y]=adams_wjl('funtest1',0,1,1,0.01,-0.5);
x
[x,y]=adams_wjl('funtest1',0,1,1,0.01,-0.5);
plot(x,y)
tic
[x,y]=adams_wjl('funtest1',0,1,1,0.01,-0.5);
time=toc
tic;[x,y]=adams_wjl('funtest1',0,1,1,0.01,-0.5);time=toc;
time
tic;[x,y]=rk_wjl('funtest1',0,1,1,0.01,-0.5);time=toc;
time
tic;[xx,yy]=rk_wjl('funtest1',0,1,1,0.01,-0.5);time=toc;
plot(x,yy-y)
hold off
plot(x,yy-y)
yyy=yy-y;
yyy
format long
yyy
format
yyy
cc
help dsolve
s= dsolve('Dy=-y +cos(2*x) - 2*sin(2*x) + 2*x*exp(-x)','y(0)=1')
normal(s)
simplify(s)
x
s1=diff(x^2*exp(-x) +cos(2*x)
s1=diff('x^2*exp(-x) +cos(2*x)')
whos
help diff
s1=diff(sym('x^2*exp(-x) +cos(2*x)'))
s
s= dsolve('Dy=-y +cos(2*t) - 2*sin(2*t) + 2*x*exp(-t)','y(0)=1')
sub(s,t,x)
syms t
s= dsolve('Dy=-y +cos(2*t) - 2*sin(2*t) + 2*t*exp(-t)','y(0)=1')
doc dsolve
s= dsolve('Dy=-y +cos(2*x) - 2*sin(2*x) + 2*x*exp(-x)','y(0)=1','x')
cc
%-- 08-5-1 上午9:31 --%
cd_temp
dir
[x,y]=adams_wjl('funtest3',0,2,1,0.2);
plot(x,y)
s= dsolve('Dy=-y +cos(2*t) - 2*sin(2*t) + 2*t*exp(-t)','y(0)=1')
s= dsolve('Dy=-y +cos(2*x) - 2*sin(2*x) + 2*x*exp(-x)','y(0)=1','x')
hold on
plot(x,cos(2*x)+exp(-x)*x^2,'*r')
plot(x,cos(2*x)+exp(-x).*x.^2,'*r')
format long
y-cos(2*x)+exp(-x).*x.^2
size(x)
size(y()
size(y)
[x,y]=adams_wjl('funtest3',0,2,1,0.2);
y-cos(2*x)+exp(-x).*x.^2
plot(x,y-cos(2*x)+exp(-x).*x.^2)
plot(x,y)
plot(x,y,'r')
clf
plot(x,y)
plot(x,y-cos(2*x)+exp(-x).*x.^2,'r')
hold on
clf
hold on
plot(x,y)
plot(x,cos(2*x)+exp(-x).*x.^2,'*r')
plot(x,y-cos(2*x)+exp(-x).*x.^2,'*y')
y
cos(2*x)+exp(-x).*x.^2
y-cos(2*x)+exp(-x).*x.^2
size(y)
size(cos(2*x)+exp(-x).*x.^2)
y-(cos(2*x)+exp(-x).*x.^2)
plot(x,y-(cos(2*x)+exp(-x).*x.^2),'^r')
[x,y]=adams_rk('funtest3',0,2,1,0.2);
[x,y]=rk_wjl('funtest3',0,2,1,0.2);
y-(cos(2*x)+exp(-x).*x.^2)
[x,y]=adams_wjl('funtest3',0,2,1,0.2);
[xx,yy]=rk_wjl('funtest3',0,2,1,0.2);
x-x
x-xx
y-yy
y-(cos(2*x)+exp(-x).*x.^2)
yy-(cos(2*x)+exp(-x).*x.^2)
printf('%20.12',yy-(cos(2*x)+exp(-x).*x.^2))
fprintf('%20.12',yy-(cos(2*x)+exp(-x).*x.^2))
fprintf('%25.15e',yy-(cos(2*x)+exp(-x).*x.^2))
x-xxclc
clc
fprintf('%25.15e\n',yy-(cos(2*x)+exp(-x).*x.^2))
fprintf('%25.23f\n',yy-(cos(2*x)+exp(-x).*x.^2))
fprintf('%+25.23f\n',yy-(cos(2*x)+exp(-x).*x.^2))
fprintf('%+25.23f\n  %+25.23f\n',y-(cos(2*x)+exp(-x).*x.^2),yy-(cos(2*x)+exp(-x).*x.^2))
fprintf('%+25.23f  %+25.23f\n',y-(cos(2*x)+exp(-x).*x.^2),yy-(cos(2*x)+exp(-x).*x.^2))
fprintf('%+25.23f  %+25.23f\n',yy-(cos(2*x)+exp(-x).*x.^2),y-(cos(2*x)+exp(-x).*x.^2))
fprintf('%+25.23f  %+25.23f\n',yy,y)
fprintf('%+25.15e  %+25.15e\n',yy,y)
plot(x,y,x,yy,'*r')
fprintf('%+25.15e ,%+25.15e\n',yy,y)
fprintf('%+25.23f,%+25.23f\n',yy,y)
y-yy
-0.36273516542313339000000-(-0.11348628197779648000000)
fprintf('%+25.23f\n',yy)
fprintf('%+25.23f\n',y)
fprintf('%+25.23f\n',cos(2*x)+exp(-x).*x.^2))
fprintf('%+25.23f\n',cos(2*x)+exp(-x).*x.^2)
[x,y]=adams_wjl('funtest3',0,2,1,0.1);
[x,y]=rk_wjl('funtest3',0,2,1,0.1);
[xx,yy]=adams_wjl('funtest3',0,2,1,0.1);
clf
plot(x,y,xx,yy,'*r')
tic;[x,y]=adams_wjl('funtest3',0,2,1,0.001);t1=toc;
t1
tic;[xx,yy]=rk_wjl('funtest3',0,2,1,0.001);t1=toc;
t1
plot(x,y,xx,yy,'*r')
cc
[x50,y50]=rk_wjl('funtest1',0,1,1,0.01,50);
[x1,y1]=rk_wjl('funtest1',0,1,1,0.01,1);
[x1f,y1f]=rk_wjl('funtest1',0,1,1,0.01,-1);
[x50f,y50f]=rk_wjl('funtest1',0,1,1,0.01,-50);
plot(x50,y50,'r-',x1,y1,'y-.',x1f,y1f,'b--',x50f,y50f,'k..')
plot(x50,y50,'-r',x1,y1,'-.y',x1f,y1f,'--b',x50f,y50f,'..k')
plot(x50,y50,'-r',x1,y1,'-.y',x1f,y1f,'--b',x50f,y50f,'..')
plot(x50,y50,'-r',x1,y1,x1f,y1f,x50f,y50f)
plot(x50,y50,x1,y1,x1f,y1f,x50f,y50f)
ylimits([-1,6e21])
ylimit([-1,6e21])
ylim([-1,6e21])
ylim([-100,6e21])
plot(x50,y50,'--',x1,y1,'-.',x1f,y1f,'-..',x50f,y50f,'-')
label
legend(['\alpha=50','\alpha=1','\alpha=-1','\alpha=-50'])
legend(['\alpha=50';'\alpha=1';'\alpha=-1';'\alpha=-50'])
legend(['\alpha=50';'\alpha=1 ';'\alpha=-1';'\alpha=-50'])
legend(['\alpha=50 ';'\alpha=1  ';'\alpha=-1 ';'\alpha=-50'])
ylim([-1,100])
50*0.01
[x2fh001,y2fh001]=rk_wjl('funtest1',0,1,1,0.01,-2);
[x2fh01,y2fh01]=rk_wjl('funtest1',0,1,1,0.1,-10);
[x2fh001,y2fh001]=rk_wjl('funtest1',0,1,1,0.01,-10);
figure
plot(x2fh01,y2fh01)
hold on
plot(x2fh001,y2fh001,'r')
[x2fh01,y2fh01]=rk_wjl('funtest1',0,1,1,0.2,-10);
plot(x2fh01,y2fh01)
[x2fh01,y2fh01]=rk_wjl('funtest1',0,1,1,0.3,-10);
plot(x2fh01,y2fh01)