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

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

[Farooq Hussain]

Thank you very much for this favour of yours.But i still have an issue with code if i have to multiply instead adding. such as given below

f1=@(x,y) A1(x)*(1.33*pi*eta*y^3)/(1+1.33*pi*eta*y^3);

kindly tell me why this vectorized formate works on for additon and subtration?

[Star Strider]

As always, my pleasure!

You will need to use element-wise operations:

f1=@(x,y) A1(x).*(1.33*pi*eta*y.^3)./(1+1.33*pi*eta*y.^3);
               ↑← HERE             ↑← AND HERE       ↑← AND HERE

That works when I run it.

I included the last two element-wise operations (multiplication and exponentiation) in my previous code. I added the element-wise multiplication by ‘A1(x)’ here, replacing the addition in the earlier version.

[Farooq Hussain]

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]

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]

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)

.