A. 11a
1a. . .
1a . . 1a .
clear all
format long
disp(' 11a')
nnp = 10; %
syms x y Dy D2y %
F = x^2+y^2+Dy^2; %
x1 = -1;
y1 = 1;
x2 = 1;
y2 = 2;
fprintf(' : F=%s\n',char(F))
fprintf(' : y(%d)=%d; y(%d)=%d\n',x1,y1,x2,y2)
dFdy = diff(F,y);
dFdy1 = diff(F,Dy);
d_dFdy1_dx = diff(dFdy1,x);
d_dFdy1_dy = diff(dFdy1,y);
d_dFdy1_dy1 = diff(dFdy1,Dy); % d(dF/dy')/dy'
dFy1dx = d_dFdy1_dx + d_dFdy1_dy*Dy + d_dFdy1_dy1*D2y;
Euler = simple(dFdy-dFy1dx);
deqEuler = [ char(Euler) '=0' ]; %
Sol = dsolve(deqEuler,'x'); %
if length(Sol)~=1 %
error(' !');
end
SolLeft = subs(Sol,x,sym(x1));
SolRight = subs(Sol,x,sym(x2));
EqLeft = [char(SolLeft) '=' char(sym(y1))];
EqRight = [char(SolRight) '=' char(sym(y2))];
Con = solve(EqLeft,EqRight); %
C1 = Con.C1;
C2 = Con.C2;
Sol1a = vpa(eval(Sol),14);
xpl = linspace(x1,x2);
y1a=subs(Sol1a,x,xpl);
11a
: F=x^2+y^2+Dy^2
: y(-1)=1; y(1)=2
2- 2- 1-
(11.4)
y ¢¢ . .
f2 = solve(deqEuler,D2y); % y''
f2 = subs (f2, {y,Dy}, {sym('y(1)'),sym('y(2)')});
rp{1} = 'function dydx = MyRightPart(x,y)';
rp{2} = 'dydx=zeros(2,1);';
rp{3} = 'dydx(1)=y(2);';
rp{4} = [ 'dydx(2)=' char(f2) ';' ];
disp(' MyRightPart.m')
fprintf ('%s\n', rp{:});
fid = fopen ('C:\Iglin\Matlab\MyRightPart.m', 'w');
fprintf (fid, '%s\n', rp{:});
fclose(fid); %
MyRightPart.m
function dydx = MyRightPart(x,y)
dydx=zeros(2,1);
dydx(1)=y(2);
dydx(2)=y(1);
(11.5)
x 1=-1; x 2=1. y 2(x 1), . y 2(x 1) : t = y 2(x 1). - , f = y 1(x 2)-2. , f t. f (t)=0. , , t . f (t) f (t)= at + b. , 2 t =0 t =1. .
|
|
xr = linspace(x1,x2,nnp+1); %
y0 = [y1,0]; % y'(x1)=0;
[xx,YY] = ode45('MyRightPart',xr,y0);
yend0 = YY(nnp+1,1)-y2
y0 = [y1,1]; % y'(x1)=1;
[xx,YY] = ode45('MyRightPart',xr,y0);
yend1 = YY(nnp+1,1)-y2
yend0 =
1.76219612254200
yend1 =
5.38905701685360
(11.3) . t = y 2(x 1) .
y0 = [y1,yend0/(yend0-yend1)]
y0 =
1.00000000000000 -0.48587364497652
. .
[xx,YY] = ode45('MyRightPart',xr,y0);
plot (xpl,y1a,'--b', xr,YY(:,1),'-r')
title ('\bfExample 11a') %
xlabel('x')
ylabel('y(x)') %
. 11.2. 11a
. .11.2 . 1a. , , , . : y ¢(x 1)=-0.48587364.
B. 11b
2. . .
2 11a. 2. .
clear all
format long
disp(' 11b')
nnp = 10; %
syms x y z Dy D2y Dz D2z %
F = Dy^2+Dz^2+2*y*z; %
x1 = -2;
y1 = 1;
z1 = 0;
x2 = 2;
y2 = 0;
z2 = 2;
fprintf(' : F=%s\n',char(F))
fprintf(' : y(%d)=%d; z(%d)=%d\n',x1,y1,x1,z1)
fprintf(' : y(%d)=%d; z(%d)=%d\n',x2,y2,x2,z2)
dFdy = diff(F,y);
dFdy1 = diff(F,Dy);
d_dFdy1_dx = diff(dFdy1,x); % d(dF/dy')/dx
d_dFdy1_dy = diff(dFdy1,y);
d_dFdy1_dy1 = diff(dFdy1,Dy);
d_dFdy1_dz = diff(dFdy1,z);
d_dFdy1_dz1 = diff(dFdy1,Dz);
dFy1dx = d_dFdy1_dx + d_dFdy1_dy*Dy + d_dFdy1_dy1*D2y + d_dFdy1_dz*Dz + d_dFdy1_dz1*D2z;
dFdz = diff(F,z);
dFdz1 = diff(F,Dz);
d_dFdz1_dx = diff(dFdz1,x); % d(dF/dz')/dx
d_dFdz1_dy = diff(dFdz1,y);
d_dFdz1_dy1 = diff(dFdz1,Dy);
d_dFdz1_dz = diff(dFdz1,z);
d_dFdz1_dz1 = diff(dFdz1,Dz);
dFz1dx = d_dFdz1_dx + d_dFdz1_dy*Dy + d_dFdz1_dy1*D2y + d_dFdz1_dz*Dz + d_dFdz1_dz1*D2z;
EulerY = simple(dFdy-dFy1dx);
EulerZ = simple(dFdz-dFz1dx);
deqEulerY = [char(EulerY) '=0']; % Y
deqEulerZ = [char(EulerZ) '=0']; % Z
Sol = dsolve(deqEulerY,deqEulerZ,'x'); %
if length(Sol)~=1 %
error(' !');
end
SolLeftY = subs(Sol.y,x,sym(x1)); % x1 y
SolLeftZ = subs(Sol.z,x,sym(x1)); % x1 z
SolRightY = subs(Sol.y,x,sym(x2)); % x2 y
SolRightZ = subs(Sol.z,x,sym(x2)); % x2 z
EqLeftY = [char(vpa(SolLeftY,14)) '=' char(sym(y1))];
EqLeftZ = [char(vpa(SolLeftZ,14)) '=' char(sym(z1))];
EqRightY = [char(vpa(SolRightY,14)) '=' char(sym(y2))];
|
|
EqRightZ = [char(vpa(SolRightZ,14)) '=' char(sym(z2))];
Con = solve(EqLeftY,EqLeftZ,EqRightY,EqRightZ);
C1 = Con.C1;
C2 = Con.C2;
C3 = Con.C3;
C4 = Con.C4;
Sol2Y = vpa(eval(Sol.y),14);
Sol2Z = vpa(eval(Sol.z),14);
xpl = linspace(x1,x2); %
y2a = subs(Sol2Y,x,xpl);
z2a = subs(Sol2Z,x,xpl);
11b
: F=Dy^2+Dz^2+2*y*z
: y(-2)=1; z(-2)=0
: y(2)=0; z(2)=2
2- 2- 4- 1- (11.1). y ¢¢, z ¢¢. y, z, y ¢, z ¢. .
f2yz = solve(deqEulerY,deqEulerZ,D2y,D2z);
f2=subs(f2yz.D2y,{y,Dy,z,Dz},{sym('y(1)'),sym('y(2)'),sym('y(3)'),sym('y(4)')});
f4=subs(f2yz.D2z,{y,Dy,z,Dz},{sym('y(1)'),sym('y(2)'),sym('y(3)'),sym('y(4)')});
rp{1} = 'function dydx = MyRightPart(x,y)';
rp{2} = 'dydx=zeros(4,1);';
rp{3} = 'dydx(1)=y(2);';
rp{4} = [ 'dydx(2)=' char(f2) ';' ];
rp{5} = 'dydx(3)=y(4);';
rp{6} = [ 'dydx(4)=' char(f4) ';' ];
disp(' MyRightPart.m')
fprintf('%s\n',rp{:});
fid = fopen ('C:\Iglin\Matlab\MyRightPart.m', 'w');
fprintf(fid,'%s\n',rp{:});
fclose(fid); %
MyRightPart.m
function dydx = MyRightPart(x,y)
dydx=zeros(4,1);
dydx(1)=y(2);
dydx(2)=y(3);
dydx(3)=y(4);
dydx(4)=y(1);
.
(11.6)
x 1 y 2(x 1)= t 1 y 4(x 1)= t 2. 2-
(11.7)
, (11.7) .
(11.8)
.
xr = linspace(x1,x2,nnp+1); %
A = zeros(2,2); %
b = zeros(2,1); %
y0 = [y1;0;z1;0]; % y'(x1)=0; z'(x1)=0;
[xx,YY] = ode45('MyRightPart',xr,y0);
b=YY(nnp+1,[1 3])' - [y2;z2]
y0 = [y1;1;z1;0]; % y'(x1)=1; z'(x1)=0;
[xx,YY] = ode45('MyRightPart',xr,y0);
A(:,1) = YY(nnp+1,[1 3])'-[y2;z2]-b;
y0 = [y1;0;z1;1]; % y'(x1)=0; z'(x1)=1;
[xx,YY] = ode45('MyRightPart',xr,y0);
A(:,2) = YY(nnp+1,[1 3])'-[y2;z2]-b
b =
13.32736606398895
11.98099902570117
A =
13.26662430099436 14.02342479230629
14.02342479230629 13.26662430099436
(11.8), . (11.1) .
yz0 = -A\b; %
y0 = [y1;yz0(1);z1;yz0(2)] %
[xx,YY] = ode45('MyRightPart',xr,y0); %
plot3(xpl,y2a,z2a,'--b',xr,YY(:,1),YY(:,3),'-r')
title ('\bfExample 11b') %
xlabel('x')
ylabel('y(x)')
zlabel('z(x)')
view(205,30)
grid on
box on
y0 =
1.00000000000000
0.42582034231673
-1.35320471612876
. 11.3. 11b
. .11.3 . 2, . : y ¢(x 1)=0.42582; z ¢(x 1)=-1.35320.
C. 11c
3. . .
(3.6) 4- 2- 2- . 4- 1-
(11.9)
t 1= y 3(x 1) t 2= y 4(x 1)
|
|
(11.10)
(3.6) , (11.10) . , , .
3 11b. 3.
clear all
format long
disp(' 11c')
nnp = 10;
syms x y Dy D2y D3y D4y %
F = D2y^2-2*Dy^2+4*y*Dy+y^2-2*y*sin(x);
x1 = -1;
y1 = 1;
Dy1 = -1;
x2 = 1;
y2 = 2;
Dy2 = 1;
fprintf(' : F=%s\n',char(F))
fprintf(' : y(%d)=%d; y''(%d)=%d\n',x1,y1,x1,Dy1)
fprintf(' : y(%d)=%d; y''(%d)=%d\n',x2,y2,x2,Dy2)
dFdy = diff(F,y);
dFdy1 = diff(F,Dy);
dFdy2 = diff(F,D2y); % dF/dy''
d_dFdy1_dx = diff(dFdy1,x); % d(dF/dy')/dx
d_dFdy1_dy = diff(dFdy1,y); % d(dF/dy')/dy
d_dFdy1_dy1 = diff(dFdy1,Dy); % d(dF/dy')/dy'
d_dFdy1_dy2 = diff(dFdy1,D2y); % d(dF/dy')/dy''
dFy1dx = d_dFdy1_dx + d_dFdy1_dy * Dy + d_dFdy1_dy1 * D2y + d_dFdy1_dy2*D3y;
d_dFdy2_dx = diff(dFdy2,x); % d(dF/dy'')/dx
d_dFdy2_dy = diff(dFdy2,y); % d(dF/dy'')/dy
d_dFdy2_dy1 = diff(dFdy2,Dy); % d(dF/dy'')/dy'
d_dFdy2_dy2 = diff(dFdy2,D2y); % d(dF/dy'')/dy''
dFy2dx = d_dFdy2_dx + d_dFdy2_dy * Dy + d_dFdy2_dy1 * D2y + d_dFdy2_dy2 * D3y;
d_dFdy2dx_dx = diff(dFy2dx,x); % d((dFy'')/dx)/dx
d_dFdy2dx_dy = diff(dFy2dx,y); % d((dFy'')/dx)/dy
d_dFdy2dx_dy1 = diff(dFy2dx,Dy); % d((dFy'')/dx)/dy'
d_dFdy2dx_dy2 = diff(dFy2dx,D2y); % d((dFy'')/dx)/dy''
d_dFdy2dx_dy3 = diff(dFy2dx,D3y); % d((dFy'')/dx)/dy'''
d2Fy2dx2 = d_dFdy2dx_dx + d_dFdy2dx_dy * Dy + d_dFdy2dx_dy1 * D2y + d_dFdy2dx_dy2 * D3y + d_dFdy2dx_dy3 * D4y;
Euler = simple(dFdy-dFy1dx+d2Fy2dx2);
deqEuler = [char(Euler) '=0']; %
Sol = dsolve (deqEuler, 'x');
if length(Sol)~=1 %
error(' !');
end
dydx = diff(Sol,x); %
slY = subs(Sol,x,sym(x1));
slDY = subs(dydx,x,sym(x1));
srY = subs(Sol,x,sym(x2));
srDY = subs(dydx,x,sym(x2));
elY = [char(vpa(slY,14)) '=' char(sym(y1))];
elDY = [char(vpa(slDY,14)) '=' char(sym(Dy1))];
erY = [char(vpa(srY,14)) '=' char(sym(y2))];
erDY = [char(vpa(srDY,14)) '=' char(sym(Dy2))];
Con = solve(elY,elDY,erY,erDY);
C1 = Con.C1;
C2 = Con.C2;
C3 = Con.C3;
C4=Con.C4;
Sol3 = vpa(eval(Sol),14); % C1-C4;
xpl = linspace(x1,x2);
y3 = subs(Sol3,x,xpl);
11c
: F=D2y^2-2*Dy^2+4*y*Dy+y^2-2*y*sin(x)
: y(-1)=1; y'(-1)=-1
: y(1)=2; y'(1)=1
yIV (11.9). , . .
f4 = solve(deqEuler,D4y); % D4y
f4 = subs(f4,{y,Dy,D2y,D3y},{sym('y(1)'),sym('y(2)'),sym('y(3)'),sym('y(4)')});
rp{1} = 'function dydx = MyRightPart(x,y)';
rp{2} = 'dydx=zeros(4,1);';
rp{3} = 'dydx(1)=y(2);';
rp{4} = 'dydx(2)=y(3);';
rp{5} = 'dydx(3)=y(4);';
rp{6} = [ 'dydx(4)=' char(f4) ';' ];
disp(' MyRightPart.m')
fprintf('%s\n',rp{:});
fid = fopen ('C:\Iglin\Matlab\MyRightPart.m', 'w');
fprintf(fid,'%s\n',rp{:});
fclose(fid); %
MyRightPart.m
function dydx = MyRightPart(x,y)
dydx=zeros(4,1);
dydx(1)=y(2);
dydx(2)=y(3);
dydx(3)=y(4);
dydx(4)=-y(1)+sin(x)-2*y(3);
(11.3, 11.10). 3 : { t 1, t 2}={0,0}; { t 1, t 2}={1,0}; { t 1, t 2}={0,1}. (11.10) . . .
|
|
xr = linspace(x1,x2,nnp+1); %
A = zeros(2,2);
b = zeros(2,1);
y0 = [y1;Dy1;0;0]; % (0,0)
[xx,YY] = ode45('MyRightPart',xr,y0); %
b = YY(nnp+1,[1 2])' - [y2;Dy2] %
y0 = [y1;Dy1;1;0]; % (1,0)
[xx,YY] = ode45('MyRightPart',xr,y0); %
A(:,1) = YY(nnp+1,[1 2])'-[y2;Dy2]-b; % 1- A
y0 = [y1;Dy1;0;1]; % (0,1)
[xx,YY] = ode45('MyRightPart',xr,y0); %
A(:,2) = YY(nnp+1,[1 2])'-[y2;Dy2]-b % 2- A
yz0 = -A\b; %
y0 = [y1;Dy1;yz0] %
[xx,YY] = ode45('MyRightPart',xr,y0); %
plot(xpl,y3,'--b',xr,YY(:,1),'-r')
title ('\bfExample 11c') %
xlabel('x')
ylabel('y(x)')
b =
-3.54609584741614
-2.69446117363805
A =
0.90929722937057 0.87079594142352
0.03850131148804 0.90929725087375
y0 =
1.00000000000000
-1.00000000000000
1.10693966258040
2.91636485466376
. 11.4. 11c
. .11.4 . 3, . : y ¢¢(x 1)=1.10694; y ¢¢¢(x 1)=2.91636.
1a, 2 3 . .