Why i don not get the graph when i add A1(x) in the function f1=@(x,y) A1(x)+(1.3​3*pi*eta*y​^3)/(1+1.3​3*pi*eta*y​^3);

2 次查看(过去 30 天)
%clc
%clear all
eta=0.25;
L=2.5;
a=0.1;
h=0.1;
b=6;
x=a:h:b;
%%%%%%%%%%%%%%%%%%
A =@(x) (0<x).*1+(0<=x & x<=L).*(1+0.5*(1-cos(2*pi*x/L)).^1/2) + (L<x).*1;
A1=@(x) (gradient(A(x))./gradient(x));
f1=@(x,y) A1(x)+(1.33*pi*eta*y^3)/(1+1.33*pi*eta*y^3);
y(1)=01;
for i=1:length(x)-1;
k1=h*f1(x(i),y(i));
k2=h*f1(x(i)+0.5*h,y(i)+.5*k1);
k3=h*f1(x(i)+.5*h,y(i)+.5*k2);
k4=h*f1(x(i)+h,y(i)+k3);
y(i+1)=y(i)+(1/6)*(k1+2*k2+2*k3+k4);
end
plot(x,y)
  4 个评论
Farooq Hussain
Farooq Hussain 2020-4-23
Thank you again very much for this latest suggestion...But when i apply the same procedure i get the same plot whether i include or remove A(x) in f3 of the code....Can you please explain
clc
clear all;
%%%%
L=1.0;
Re=1000;
k=1.4;
we=137;
d=0.8; %%%% d=sigma
alphas=0;
%%% Step size%%%
a1=0;
h=0.05;
n=5;
x=[a1:h:n];
%%% Functions %%%
A = @(x)(0<x).*1+(0<=x & x<=L).*(1+0.5*(1-cos(2*pi*x/L))).^1/2 + (L<x).*1;
f1=@(x,u,v,cp,H,R)v;
f2=@(x,u,v,cp,H,R)H;
f3=@(x,u,v,cp,H,R)(-1+alphas+alphas*R.^3)+A(x).*u;
f4=@(x,u,v,cp,H,R)(2*u*v*(1-alphas))/(-1+alphas+alphas*R^3);
f5=@(x,u,v,cp,H,R)(2*(R^(-3*k)-(R^-1))/we+0.5*d*(-1+R^(-3*k))-0.5*cp-u*v*R*H-1.5*(u^2)*(H^2)-(4*u*R*H)/(Re))./(R*u^2);
%%% Initial Conditions %%%
cp(1)=0;
u(1)=01;
v(1)=0.0154;
R(1)=01;
H(1)=0.021;
for i=1:length(x)-1;
m1=h*f1(x(i), u(i), v(i), cp(i), H(i), R(i));
n1=h*f2(x(i), u(i), v(i), cp(i), H(i), R(i));
p1=h*f3(x(i), u(i), v(i), cp(i), H(i), R(i));
q1=h*f4(x(i), u(i), v(i), cp(i), H(i), R(i));
r1=h*f5(x(i), u(i), v(i), cp(i), H(i), R(i));
m2=h*f1(x(i)+0.5*h, u(i)+0.5*m1, v(i)+0.5*n1, cp(i)+0.5*p1, H(i)+0.5*q1, R(i)+0.5*r1);
n2=h*f2(x(i)+0.5*h, u(i)+0.5*m1, v(i)+0.5*n1, cp(i)+0.5*p1, H(i)+0.5*q1, R(i)+0.5*r1);
p2=h*f3(x(i)+0.5*h, u(i)+0.5*m1, v(i)+0.5*n1, cp(i)+0.5*p1, H(i)+0.5*q1, R(i)+0.5*r1);
q2=h*f4(x(i)+0.5*h, u(i)+0.5*m1, v(i)+0.5*n1, cp(i)+0.5*p1, H(i)+0.5*q1, R(i)+0.5*r1);
r2=h*f5(x(i)+0.5*h, u(i)+0.5*m1, v(i)+0.5*n1, cp(i)+0.5*p1, H(i)+0.5*q1, R(i)+0.5*r1);
m3=h*f1(x(i)+0.5*h, u(i)+0.5*m2, v(i)+0.5*n2, cp(i)+0.5*p2, H(i)+0.5*q2, R(i)+0.5*r2);
n3=h*f2(x(i)+0.5*h, u(i)+0.5*m2, v(i)+0.5*n2, cp(i)+0.5*p2, H(i)+0.5*q2, R(i)+0.5*r2);
p3=h*f3(x(i)+0.5*h, u(i)+0.5*m2, v(i)+0.5*n2, cp(i)+0.5*p2, H(i)+0.5*q2, R(i)+0.5*r2);
q3=h*f4(x(i)+0.5*h, u(i)+0.5*m2, v(i)+0.5*n2, cp(i)+0.5*p2, H(i)+0.5*q2, R(i)+0.5*r2);
r3=h*f5(x(i)+0.5*h, u(i)+0.5*m2, v(i)+0.5*n2, cp(i)+0.5*p2, H(i)+0.5*q2, R(i)+0.5*r2);
m4=h*f1(x(i)+h, u(i)+m3, v(i)+n3, cp(i)+p3, H(i)+q3, R(i)+r3);
n4=h*f2(x(i)+h, u(i)+m3, v(i)+n3, cp(i)+p3, H(i)+q3, R(i)+q3);
p4=h*f3(x(i)+h, u(i)+m3, v(i)+n3, cp(i)+p3, H(i)+q3, R(i)+r3);
q4=h*f4(x(i)+h, u(i)+m3, v(i)+n3, cp(i)+p3, H(i)+q3, R(i)+r3);
r4=h*f5(x(i)+h, u(i)+m3, v(i)+n3, cp(i)+p3, H(i)+q3, R(i)+r3);
u(i+1)=u(i)+(1/6)*(m1+2*m2+2*m3+m4);
v(i+1)=v(i)+(1/6)*(n1+2*n2+2*n3+n4);
cp(i+1)=cp(i)+(1/6)*(p1+2*p2+2*p3+p4);
H(i+1)=H(i)+(1/6)*(q1+2*q2+2*q3+q4);
R(i+1)=R(i)+(1/6)*(r1+2*r2+2*r3+r4);
end
plot(x,u,'k')
Star Strider
Star Strider 2020-4-23
The ‘f3’ function likely makes a very small contribution, regardless. (The code is so complicated that it is difficult to follow.)
Plotting ‘f3’ alone (after the end of the loop so that all the vectors are present) with:
f3=@(x,u,v,cp,H,R)(-1+alphas+alphas*R.^3)+A(x).*u;
or without:
f3=@(x,u,v,cp,H,R)(-1+alphas+alphas*R.^3);%+A(x).*u;
the ‘A(x)’ term definitely makes a difference:
f3v = f3(x+h, u+m3, v+n3, cp+p3, H+q3, R+r3);
figure
plot(x, h*f3v);
So ‘f3’ works correctly.
.

请先登录,再进行评论。

采纳的回答

Star Strider
Star Strider 2020-4-18
Do not use the for loop here. Let vectorisation work for you.
Try this:
eta=0.25;
L=2.5;
a=0.1;
h=0.1;
b=6;
x=a:h:b;
%%%%%%%%%%%%%%%%%%
A =@(x) (0<x).*1+(0<=x & x<=L).*(1+0.5*(1-cos(2*pi*x/L)).^1/2) + (L<x).*1;
A1=@(x) (gradient(A(x))./gradient(x));
f1=@(x,y) A1(x)+(1.33*pi*eta*y.^3)./(1+1.33*pi*eta*y.^3);
y(1)=01;
% for i=1:length(x)-1;
k1=h*f1(x,y);
k2=h*f1(x+0.5*h,y+.5*k1);
k3=h*f1(x+.5*h,y+.5*k2);
k4=h*f1(x+h,y+k3);
y = +(1/6)*(k1+2*k2+2*k3+k4);
% end
plot(x,y)
.

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Directed Graphs 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by