You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
Please, I would appreciate any assistance on debugging the attached code.
2 views (last 30 days)
Show older comments
The first section of the code is okay, but I am having it difficult solving for uu and ud in the attached code. I want to continue solving for the values of uu and ud as long the conditions in the while loops are valid. Pls, see the attached file.
This is the pseudo code:
for j=1:N
k=1;
uu(k,j)=0;
Lou(k,j)=2;
p(k,j)=6000;
pQ(k,j)=asin(1+Lou(k,j));
Lh2(j)=200;
k22(k,j)=Lh2(j)./Lou(k,j);
if p(k,j)>pQ(k,j)
while Lou(k,j)<Lh2
k22f(k,j)=Lh2(j)./(Lou(k,j)+uus);
funLu(k,j)=asin(1+k22f(k,j));
testfunLu=@(x)eval(subs(funLu(k,j),uus,x))-p(k,j);
uuf(k,j)=fsolve(testfunLu,55);
uu(k+1,j)=uuf(k,j);
Lou(k+1,j)=Lou(k,j)+uu(k+1,j);
k22(k+1,j)=Lh2(j)./(Lou(k+1,j)+uu(k+1,j));
p(k+1,j)=p(k,j);
k=k+1
end
end
end
9 Comments
Isaac
on 1 Apr 2016
Edited: Walter Roberson
on 2 Apr 2016
This is the part of the code, that has bug:
if p(k,j)>=pQu1(k,j)| p(k,j)>=pQd1(k,j) %displacements within the host formation (Case 1)
while Lou(k,j)<Lh2 & Lod(k,j)<Lh2
%Numerically Solve for ud at t=k+1 using Newton Raphson
funLu(k,j)=KICfm(j)./sqrt(pi*(Lou(k,j)+uus))+S2+0.5*grav*(Lou(k,j)+uus)*rhom; %up-tip
uuf{k,j}=real(double(solve(funLu(k,j)==p(k,j),uus)));
uu(k+1,j)=uuf{k,j}(2);
if uu(k+1,j)<0
uu(k+1,j)=0;
end
funLd(k,j)=KICfm(j)./sqrt(pi*(Lod(k,j)+uds))+S2-0.5*grav*(Lod(k,j)+uds)*rhom; %down-tip
udf{k,j}=real(double(solve(funLd(k,j)==p(k,j),uds)));
ud(k+1,j)=udf{k,j}(1);
if ud(k+1,j)<0
ud(k+1,j)=0;
end
Lou(k+1,j)=Lou(k,j)+uu(k+1,j);
pQu1(k+1,j)=KICfm(j)./sqrt(pi*Lou(k+1,j))+S2+0.5*grav*Lou(k+1,j)*rhom;
Lod(k+1,j)=Lod(k,j)+ud(k+1,j);
pQd1(k+1,j)=KICfm(j)./sqrt(pi*Lod(k+1,j))+S2-0.5*grav*Lod(k+1,j)*rhom;
p(k+1,j)=p(k,j);
%Case 2
pQd2(k+1,j)= pQd2(k,j);
pQu3(k+1,j)=pQu3(k,j);
k=k+1;
end
if p(k,j)>pQd2(k,j)& Lod(k,j)>Lh2 & Lou(k,j)<Lh2 %case-2: up-tip within host formation and down-tip outside
Lod(k,j)=Lh2(1);
%Lou(k,j)=Lh2(1);
while (Lou(k,j)<=Lh2 & Lh2<Lod(k,j)<Lh3 )
%Calculating the displacements
%lower tip
%down-tip critical pressure function at Q
k22f(k,j)=Lh2(j)./(Lod(k,j)+uds);
funLd(k,j)=(sqrt(pi)*KICd(j)./(Lod(k,j)+uds)-grav*(Lod(k,j)+uds)*pi*rhom-sqrt((1+k22f(k,j))/(1-k22f(k,j)))*S2+k22f(k,j).*sqrt((1+k22f(k,j))/(1-k22f(k,j)))*S2+sqrt(1-k22f(k,j).^2)*S3+pi*S3+grav*(Lod(k,j)+uds)*rhom*asin(sqrt((1-k22f(k,j))/2))+(grav*(Lod(k,j)+uds)*rhom+2*S2-2*S3).*asin(sqrt((1+k22f(k,j))/2)))./(-sqrt((1+k22f(k,j))/(1-k22f(k,j)))+k22f(k,j).*sqrt((1+k22f(k,j))/(1-k22f(k,j)))+sqrt(1-k22f(k,j).^2)+pi);%down-tip
udf{k,j}=real(double(solve(funLd(k,j)==p(k,j),uds)));
ud(k+1,j)=udf{k,j}(2);
if ud(k+1,j)<0
ud(k+1,j)=0;
end
%upper tip
k22(k,j)=Lh2(j)./(Lo(k,j)+ud(k,j));
funLu(k,j)=-(2*KICfm(j)*sqrt(pi)+sqrt((Lou(k,j)+uus)).*(grav*(Lou(k,j)+uus)*rhom*pi+2*(sqrt(1-k22(k,j).^2)*S2-(sqrt((1-k22(k,j))/(1+k22(k,j)))+k22(k,j).*sqrt((1-k22(k,j))/(1+k22(k,j)))-pi)*S3))+4*sqrt((Lou(k,j)+uus))*(S2-S3).*asin(sqrt((1+k22(k,j))/2)))./(2*sqrt((Lou(k,j)+uus)).*(sqrt((1-k22(k,j))/(1+k22(k,j)))+k22(k,j).*sqrt((1-k22(k,j))/(1+k22(k,j)))-sqrt(1-k22(k,j).^2)-pi)); %up-tip
uuf{k,j}=real(double(solve(funLu(k,j)==p(k,j),uus)));
uu(k+1,j)=uuf{k,j}(2);
if uu(k+1,j)<0
uu(k+1,j)=0;
end
Lod(k+1,j)=Lod(k,j)+ud(k+1,j);
k22(k+1,j)=Lh2(j)./(Lod(k+1,j)+ud(k+1,j));
pQd2(k+1,j)=(sqrt(pi)*KICd(j)./Lod(k+1,j)-grav*Lod(k+1,j)*pi*rhom-sqrt((1+k22(k+1,j))/(1-k22(k+1,j)))*S2+k22(k+1,j).*sqrt((1+k22(k+1,j))/(1-k22(k+1,j)))*S2+sqrt(1-k22(k+1,j).^2)*S3+pi*S3+grav*Lod(k+1,j)*rhom*asin(sqrt((1-k22(k+1,j))/2))+(grav*Lod(k+1,j)*rhom+2*S2-2*S3).*asin(sqrt((1+k22(k+1,j))/2)))./(-sqrt((1+k22(k+1,j))/(1-k22(k+1,j)))+k22(k+1,j).*sqrt((1+k22(k+1,j))/(1-k22(k+1,j)))+sqrt(1-k22(k+1,j).^2)+pi);%down-tip
Lou(k+1,j)=Lou(k,j)+uu(k+1,j);
pQu2(k+1,j)=-(2*KICfm(j)*sqrt(pi)+sqrt(Lou(k+1,j)).*(grav*Lou(k+1,j)*rhom*pi+2*(sqrt(1-k22(k+1,j).^2)*S2-(sqrt((1-k22(k+1,j))/(1+k22(k+1,j)))+k22(k+1,j).*sqrt((1-k22(k+1,j))/(1+k22(k+1,j)))-pi)*S3))+4*sqrt(Lou(k+1,j))*(S2-S3).*asin(sqrt((1+k22(k+1,j))/2)))./(2*sqrt(Lou(k+1,j)).*(sqrt((1-k22(k+1,j))/(1+k22(k+1,j)))+k22(k+1,j).*sqrt((1-k22(k+1,j))/(1+k22(k+1,j)))-sqrt(1-k22(k+1,j).^2)-pi)); %up-tip
p(k+1,j)=p(k,j);
%Case 3
pQu3(k+1,j)= pQu3(k,j);
k=k+1;
end
elseif p(k,j)>pQu3(k,j)& Lou(k,j)>Lh2 & Lod(k,j)<Lh2 %case-3: down-tip within host formation and up-tip outside
%Lod(k,j)=Lh2(1);
Lou(k,j)=Lh2(1);
while(Lod(k,j)<=Lh2 & Lh2<Lou(k,j)<Lh3)
%Calculating the displacements
%upper tip
k13f(k,j)=Lh2(j)./(Lou(k,j)+uus);
funLu(k,j)=(2*KICu(j)*sqrt(pi)+sqrt((Lou(k,j)+uus)).*(grav*(Lou(k,j)+uus)*rhom*pi+2*(sqrt(1-k13f(k,j).^2)*S1+(-sqrt((1+k13f(k,j))./(1-k13f(k,j)))+k13f(k,j).*sqrt((1+k13f(k,j))./(1-k13f(k,j)))+pi)*S2))+4*sqrt((Lou(k,j)+uus))*(S1-S2).*asin(sqrt((1-k13f(k,j))/2)))./(2*sqrt((Lou(k,j)+uus)).*(-sqrt((1+k13f(k,j))./(1-k13f(k,j)))+k13f(k,j).*sqrt((1+k13f(k,j))./(1-k13f(k,j)))+sqrt(1-k13f(k,j).^2)+pi)); %up-tip
testfunLu=@(x)eval(subs(funLu(k,j),uus,x))-p(k,j);
%uuf{k,j}=real(double(fsolve(funLu(k,j)==p(k,j),uus)));
uuf(k,j)=fsolve(testfunLu,55);
uu(k+1,j)=uuf(k,f);
if uu(k+1,j)<0
uu(k+1,j)=0;
end
%lower tip
k13(k,j)=Lh2(j)./(Lou(k,j)+uu(k,j));
funLd(k,j)=0.5./(sqrt(Lod(k,j))*pi).*(2*KICfm(j)*sqrt(pi)+sqrt((Lod(k,j)+uds)).*(-grav*(Lod(k,j)+uds)*pi*rhom+2*(-sqrt((1-k13(k,j))/(1+k13(k,j))).*(1+k13(k,j))*S1+(sqrt((1-k13(k,j))/(1+k13(k,j)))+k13(k,j).*sqrt((1-k13(k,j))/(1+k13(k,j)))+pi)*S2))+4*sqrt((Lod(k,j)+uds))*(S1-S2).*asin(sqrt((1-k13(k,j))/2))); %down-tip
testfunLd=@(x)eval(subs(funLd(k,j),uds,x))-p(k,j);
udf(k,j)=fsolve(testfunLd,55);
ud(k+1,j)=udf(k,f);
if ud(k+1,j)<0
ud(k+1,j)=0;
end
Lou(k+1,j)=Lou(k,j)+uu(k+1,j);
k13(k+1,j)=Lh2(j)./(Lou(k+1,j)+uu(k+1,j));
pQu3(k+1,j)=(2*KICu(j)*sqrt(pi)+sqrt((Lou(k+1,j))).*(grav*(Lou(k+1,j))*rhom*pi+2*(sqrt(1-k13(k+1,j).^2)*S1+(-sqrt((1+k13(k+1,j))/(1-k13(k+1,j)))+k13(k+1,j).*sqrt((1+k13(k+1,j))/(1-k13(k+1,j)))+pi)*S2))+4*sqrt((Lou(k+1,j)))*(S1-S2).*asin(sqrt((1-k13(k+1,j))/2)))./(2*sqrt((Lou(k+1,j))).*(-sqrt((1+k13(k+1,j))/(1-k13(k+1,j)))+k13(k+1,j).*sqrt((1+k13(k+1,j))/(1-k13(k+1,j)))+sqrt(1-k13(k+1,j).^2)+pi)); %up-tip
Lod(k+1,j)=Lod(k,j)+ud(k+1,j);
pQd3(k+1,j)=0.5./(sqrt(Lod(k+1,j))*pi).*(2*KICfm(j)*sqrt(pi)+sqrt(Lod(k+1,j)).*(-grav*Lod(k+1,j)*pi*rhom+2*(-sqrt((1-k13(k+1,j))/(1+k13(k+1,j))).*(1+k13(k+1,j))*S1+(sqrt((1-k13(k+1,j))/(1+k13(k+1,j)))+k13(k+1,j).*sqrt((1-k13(k+1,j))/(1+k13(k+1,j)))+pi)*S2))+4*sqrt(Lod(k+1,j))*(S1-S2).*asin(sqrt((1-k13(k+1,j))/2))); %down-tip
p(k+1,j)=p(k,j);
%Case 4
pQu4(k+1,j)= pQu4(k,j);
pQd4(k+1,j)= pQd4(k,j);
k=k+1;
end
end
end
end
Chad Greene
on 1 Apr 2016
Yikes. I can tell right away that you're unlikely to get a response because there's too much to sift through. Please trim your code down to a MWE and clearly describe what you mean when you say you're having difficulty.
Isaac
on 1 Apr 2016
Edited: per isakson
on 2 Apr 2016
This is a summary:
for j=1:N
k=1
uu(k,j)=0 %initialization
Lou(k,j)=2; %initialization
while Lou(k,j)<L
uu(k+1,j)=fsolve...
.
.
The problem is how to use fsolve in solving for the value of uu for each time it goes into the while loop; k is the counter
jgg
on 1 Apr 2016
This isn't a minimal example. How is Lou(k,j) updated? What is the solving criterion? You're making this very difficult for anyone to help you with.
For instance, the code you've provided doesn't even run because uuf{k,j}(2) is a single value, not a vector. Additionally, you're not pre-declaring any of your variables, so they change size every iteration which makes it impossible to tell how anything should be structured.
Isaac
on 1 Apr 2016
k is a counter for the while loop. I am solving for uu as long as the condition in the while loop is valid. So the length of k depends on the iteration steps.
Lou is updated in the while loop.
Lou(k+1,j)=Lou(k,j)+uu(k,j)
My
Isaac
on 1 Apr 2016
Edited: per isakson
on 2 Apr 2016
I have tried using the fsolve this way also
testfunLu=@(x)eval(subs(funLu(k,j),uus,x))-p(k,j);
uuf(k,j)=fsolve(testfunLu,55);
Isaac
on 1 Apr 2016
Edited: per isakson
on 2 Apr 2016
This is the pseudo code:
for j=1:N
k=1;
uu(k,j)=0;
Lou(k,j)=2;
p(k,j)=6000;
pQ(k,j)=asin(1+Lou(k,j));
Lh2(j)=200;
k22(k,j)=Lh2(j)./Lou(k,j)
if p(k,j)>pQ(k,j)
while Lou(k,j)<Lh2
k22f(k,j)=Lh2(j)./(Lou(k,j)+uus);
funLu(k,j)=asin(1+k22f(k,j));
testfunLu=@(x)eval(subs(funLu(k,j),uus,x))-p(k,j);
uuf(k,j)=fsolve(testfunLu,55);
uu(k+1,j)=uuf(k,j);
Lou(k+1,j)=Lou(k,j)+uu(k+1,j);
k22(k+1,j)=Lh2(j)./(Lou(k+1,j)+uu(k+1,j));
p(k+1,j)=p(k,j);
k=k+1
end
end
end
Answers (0)
See Also
Categories
Find more on Loops and Conditional Statements in Help Center and File Exchange
Tags
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)