Why i don not get the graph when i add A1(x) in the function f1=@(x,y) A1(x)+(1.33pietay^3)/(1+1.33pietay^3);

Question

Farooq Hussain 2020-4-18

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/518918-why-i-don-not-get-the-graph-when-i-add-a1-x-in-the-function-f1-x-y-a1-x-1-33-pi-eta-y-3-1-1

评论： Star Strider 2020-4-23

%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 个评论
显示 2更早的评论隐藏 2更早的评论

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 2020-4-23

在 MATLAB Online 中打开

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.

.

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2020-4-18

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/518918-why-i-don-not-get-the-graph-when-i-add-a1-x-in-the-function-f1-x-y-a1-x-1-33-pi-eta-y-3-1-1#answer_426820

在 MATLAB Online 中打开

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 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Why i don not get the graph when i add A1(x) in the function f1=@(x,y) A1(x)+(1.33pietay^3)/(1+1.33pietay^3);

4 个评论
显示 2更早的评论隐藏 2更早的评论

采纳的回答

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

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);

4 个评论 显示 2更早的评论隐藏 2更早的评论

采纳的回答

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

Why i don not get the graph when i add A1(x) in the function f1=@(x,y) A1(x)+(1.33pietay^3)/(1+1.33pietay^3);

4 个评论
显示 2更早的评论隐藏 2更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论