Matlab: Solving Boundary Value Problems

[@brynntheoose]

This was extremely helpful, especially the section on defining the boundary conditions as a vector, thank you!

[@sreeeeee777]

THANK YOU SO MUCH. It was very wonderfully explained. Felt so easy about BVPs. Hats off :)

[@tasslebear82]

Hi jake, thanks alot for your prompt response. my equation is actually linear in y, its just the coefficients of y, y' and y" are dependant on my independant variable u. i've since looked at, and printed off (paper is so much easier on the eyes) that tutorial by Kierzenka, thanks alot for suggesting it. they do quite a few examples that mirror my own problem to an extent.. i imagine i'll have the code working pretty soon. thanks again!!

[@azizi0407]

Thank you for this tutorial. It does help me a lot!

[@brynntheoose]

Just curious, how does changing the initial mesh to having just 5 points affect your final answer?

[@edholanda5897]

Professor how would I solve this system of first order edos numerically by plotting the graph for the different values ​​of (n). the derivatives are in relation to ha (r). a'/r = -e^2*v^2*(g^2 - 1) g' = - a*g/r given the boundary conditions a(0) = n a(inf)=0 g(0) = 0 g(inf)=1 where (e)=0.5 and (v)=1 are constant and n variavel. please give a helping hand there, I looked for and did not find any problems like this on the matlab website.

[@jakeblanchard]

I think you want to treat this like an initial value problem and guess values for n such that you get the right values at t=infinity. One concern is that treating n as a variable might help you get the right value for either a(inf) or g(inf), but I'm skeptical you can get both. Play around with that approach and see if it gets you anywhere. If the system behaves itself, you might be able to get a value for n by treating the whole thing as a root-finding problem.

[@anishdasgupta100]

I was just wondering if there is a built in function for Fourier constants. Otherwise I will have to define the function as an M-file with n as a parameter and integrate it numerically, right?But I usually solve for Fourier constants using orthogonality condition so there will 2 integration per An. Also what is the Matlab func for finding eigen values? Do you have any tutorials regarding this? Thanking You.

[@jakeblanchard]

I've never solved a problem like that. I would need more details to be of any help.

[@uoweme5grand]

I got an error saying: ??? Error using ==> bvp4c Unable to solve the collocation equations -- a singular Jacobian encountered. how do i solve this?

[@marionaranjo76]

Thanks man, this video was super helpful!!

[@welidbenchouche]

very awesome, however, i find Unable to solve the collocation equations -- a singular Jacobian encountered. how to overcome this error, for the record, my ode is a first order eq

[@jakeblanchard]

There are a variety of things that can cause a singular Jacobian. You might try a different initial guess before solving. The following link has a thorough discussion of the possible causes for this error. www.mathworks.com/matlabcentral/answers/120468-why-do-i-receive-an-error-about-a-singular-jacobian-when-i-use-the-bvp4c-function-within-matlab

[@jakeblanchard]

By the way, do you really have a single, first-order equation? I don't think that will work. You would need at least two first order equations for a boundary value problem. Feel free to send me a description of the problem so that I can comment in more detail.

[@ciwstevie]

Very nice explanations on solving boundary value problems! THANKS FOR THAT! But when trying to compile the code in Matlab (R2008a) I realized that there are semicolons (;) missing for the vectors dydx and res in the functions describing the differential equation and the boundary conditions, respectively. Matlab needs vertical vectors to handle the ODE. Or am I wrong? At least it worked out fine when adding the semicolons.

[@aliexpress.official]

thank you so much. that function is the most needlessly convoluted rubbish

[@abdullahokyar8509]

Hi,When I write on MATLAB ,I I see that error ??? Undefined function or method 'linspace' for input arguments of type 'double'. Error in ==> bvp5 at 3 solinit=bvpinit(linspace(xlow,xhigh,10),(1-1));

[@jakeblanchard]

That seems strange. Maybe you could send me your script. [e-mail]

[@abdullahokyar8509]

I send it

[@anishdasgupta100]

Dear Sir, I want to solve a PDE using eigen values and separation of variables. The final solution is a Fourier series. Actually, it is a 2 directional flow of fluid problem. can you plz help me with that?I can solve the problem but am not sure as to the best way to calculate the Fourier constants (An). Also is there a way to find the eigenvalue (as a function of n) using matlab?

[@tasslebear82]

hi, this is a really excellent video ya put up, very well explained!! im using matlab to solve a 2nd order BVP at the mo, however my system of first order ode's have non constant coefficients dependant on my independant variable u. is it just a matter of defining the coeffs and multiplying them by my y(1) and y(2) in my ode file before i execute the script? i seem to be getting errors relating back to my ode file. my bc's consist of y and y' at one given point. any insight wud be great, thanks

[@sonnyskaa]

In regards to picking the initial guess you pick y=1 and y'=-1. What is the reasons for these values? Since bvp4c uses collocation points do we need to satisfy the boundary conditions with our choice of y and y'?

[@jakeblanchard]

In general, I would suggest that for any system you are solving, you should understand it well enough, prior to solving numerically, that you have a good idea what y and y' might be. In this case, I knew that y(0)=1 so I guessed that all values of y will be of order 1 so I used that as a guess. For y', there was nothing in the differential equation to lead me to expect a large gradient, so order 1 for that was also sufficient. I guessed that it would be negative since the boundary value on the right was lower than the one on the left.

[@sonnyskaa]

Thanks for the help. I know quite a bit about the solution. I know that the function is monotonically decreasing (hence picking a negative value for y'). I have a bound on y' (between 0 and a a minus constant depending on parameters). Since the left BC is y(0)=0 and the function is monotonically decreasing I picked y=0 and y'=-1. This seemed reasonable. After looking into what people picked as initial guess people stated that they picked guess for y that satisfies both boundary conditions and then differentiated it to have guess for y'. Thus I was a little confused, but it seemed that my original approach (picking y=0 and y'=-1) should work. Thanks for you help

[@jakeblanchard]

Sonny Skaaning Some problems require better guesses than others. For linear problems it generally isn't critical.

[@SuperExodus13]

Just a question to make sure I understand this correctly. If I have a third order nonlinear ode do I need to make new variables for each of my derivatives? Like how you say z=dy/dx, but I would also need lets say a=d2y/dx2, b=d3y/dx3? I have 1 Neuman BC, and 2 Diriclet BC's.

[@SuperExodus13]

+Nathaniel Werner Also, my range is from -.09 to 2

[@jakeblanchard]

+Nathaniel Werner For a third order equation, you want to reduce it to three first order equations. So your variables could be y, a, and b, with a=dy/dx and b=d2y/dx2.

[@deadrunner983]

Thanks a million!

[@hearttoolove]

yes, I have been getting the same problem typing the exact script....there must be some problem I feel

[@sibyjoseplathottam4828]

Very helpful video.

[@maahir4763]

Thanks a lot great video.

[@tslee14]

Thanks so much for this!

[@safdar9688]

help me sir please this cod cannot run (Error in newpaperssss (line 18) sol=bvp4c(@bvpexam2,@bcexam2,solinit); )

[@Fridaydude]

Good day sir. In regards to picking appropriate boundary conditions for a problem, do you always have to consider boundary conditions that are functions of lower rank ( derivative-wise ) in comparison with the problem itself? I.e, for your first example in the video, which is a 2nd order derivative equation, do boundary conditions have to include 1st derivative and zero order terms ( namely y' and y )? Or you can use y'' as well as boundary conditions? Thanks in advance.

[@jakeblanchard]

+Paraskevas Dimitrios The short answer is: yes. Every problem I've ever done was as you described. But, I suppose there might be situations which are less straightforward.

[@Fridaydude]

+Jake Blanchard Thanks a lot, appreciated

[@thangdovuminh]

thank you. But if i want to receive the solution of z, how should i do?

[@jakeblanchard]

plot(xint, Sxint(2,:))

[@thangdovuminh]

@@jakeblanchard thanks

[@thangdovuminh]

is there any rules to choose the initial guess point? or just pick random ?

[@jakeblanchard]

@@thangdovuminh No, there really aren't. I usually try to focus on the scientific problem you are trying to solve and estimate a solution that way. Sometimes you can just focus on the equations, simplify them in a way that doesn't affect the solution much, and then use the solution to the simplified solution as an initial guess for the whole solution. Mostly, you just need to be creative.

[@thangdovuminh]

@@jakeblanchard thank you so much

[@Abrecadebre]

Thank you so much for the tutorial. However, when I'm trying to solve this (d2y/dx2 + (1/x)(dy/dx) = -100-100x^2), this warning always pops out when i set my xlow=0 Error using bvp4c (line 251) Unable to solve the collocation equations -- a singular Jacobian encountered. Error in bvp5 (line 5) sol=bvp4c(@bvp5ode,@bvp5bc,solinit); how can i solve this? Thanks.

[@jakeblanchard]

You might try bvp5c, since it's supposed to do better with singular equations. Also, take a look at the discussion of Singular BVPs here: www.mathworks.com/help/matlab/ref/bvpset.html

[@vamsikrishna7194]

Hey how can i set boundary conditions(w.r.t space only) along with the initial conditions(w.r.t. time only) for an O.D.E?

[@jakeblanchard]

If you want to include both initial conditions and initial conditions, then it sounds to me like you have a PDE, not an ODE. Depending on the equation, pdepe might be valuable. This is a Matlab routine for parabolic and elliptic equations.

[@vamsikrishna7194]

My problem actually is a composite wall system whose end temperatures are known and the whole system has a uniform temperature at t=0 ; and to get a temperature vs. distance plot. So I cant rely on pdepe as I discretized the space using finite difference method!

[@jakeblanchard]

So it sounds like you just want the steady state solution. Then you ignore the initial conditions and solve the boundary value problem.

[@vamsikrishna7194]

It's a time derivative and I have to plot for like various times...so it's not a steady state problem, instead it is an unsteady state problem with 2 boundary conditions given(at the ends of the composite wall)

[@jakeblanchard]

OK, then it isn't an ODE and you can't use bvp4c. I would use pdepe, but you say that you've already discretized things. I would still try pdepe (it will difference/discretize for you), but that is up to you. If you have already differenced things, then you just need to factor in the boundary and initial conditions and solve things. The solution approach will depend on whether you used implicit or explicit differencing.

[@hammaminadjima3813]

Thanks a lot ...i ask how can i solve this problem : MATLAB : Jacobi : singoler ?

[@jakeblanchard]

Any number of things could cause this error. It could be a bad initial guess, or an ill-posed problem, or something else. I would have to see your code and a description of the problem in order to say any more.

[@hammaminadjima3813]

@@jakeblanchard Ok .Can you give me your email, i will email you

[@jakeblanchard]

@@hammaminadjima3813 blanchard "at" engr.wisc.edu

[@sonnyskaa]

If you change the right BC to y(pi)=0 hence changing xhigh=pi. When I change the mesh in bvpinit between 5,10,100 etc. I get the same graph, but EXTREMELY different values for the y-coordinate. Why is that??? You Sxint = deval(sol,xint) gives a 2:mess vector. Why is it 2D and what does the values indicate (thus what does the values in row 1 and row 2 mean) Thanks

[@sonnyskaa]

I believe the Sxint gives a 2D vectore with solutions for y(1) and y(2) if I am not mistaken. I still dont understand why I get EXTREMELY different values for y(1) and y(2) when I change Right BC to pi and change the mesh in bvpinit.

[@jakeblanchard]

Yes, that new boundary condition changes everything. Try this: solve it with y(1.01 pi)=0. Then try y(1.001 pi)=0. See what happens as you get closer to y(pi)=0.

[@jakeblanchard]

Sonny Skaaning Regarding the 2-D nature of the problem, the way I set up the solution one column holds y and the other holds y'

[@sonnyskaa]

Jake Blanchard Yes I figured that out after some inventigation. But still dont understand that if I change the mess (in linspace(xlow,xhigh,mess) then I get different y,y' results with Right BC set at pi. I have been thinking about our initial guess vector. Is this guess vector a guess at each mess point or just the Left BC?

[@jakeblanchard]

Sonny Skaaning You guess everything at every mesh point. But the issue with the mesh spacing goes back to the solution of the equation. My math is rusty, but it appears to me that the solution to the problem you proposed is singular. Hence, the finer the mesh the finer you resolve the singularity. You should study that problem analytically before you try it numerically.

[@baacademy1166]

Thank you very much for this wonderful tutorial. Please help me with this question: When one of the boundary conditions is dydx=0. why 1-y(1) become -y(1)? Thank you

[@jakeblanchard]

I'm not sure I understand your question, so feel free to ask again if my answer doesn't help. The example in this video does not have any boundary conditions involving dydx. If you wanted dydx=0 at x=a, you would need something like ya(2) as a term in the res term of the bvp4bc function.

[@baacademy1166]

@@jakeblanchard thank you for the response sir. it became clear watching the video again. but receive the message: Error using bvp4c Unable to solve the collocation equations -- a singular Jacobian encountered .when I tried running my code below. plz help note: it run well when I exclude the term -(2/r)*Cs(2) function bvp4 rlow=0; rhigh=1; B=1; T=2; solinit=bvpinit(linspace(rlow,rhigh,10),[1 1]); sol=bvp4c(@bvp4ode, @bvp4bc, solinit); rint=linspace(rlow,rhigh); Srint=deval(sol,rint); plot(rint,Srint(1,:)) xlabel('radius') ylabel('concentration') title('substrate profile when reaction dominate') function dCsdr=bvp4ode(r,Cs) dCsdr=[Cs(2) -(2/r)*Cs(2)+((9*T^2) *(Cs(1)/(1+B*Cs(1))))]; end end function res=bvp4bc(Csa,Csb) res=[Csa(1)-1 Csb(2)]; end

[@jakeblanchard]

@@baacademy1166 Two things. You have to define B and T inside the bvp4ode function. More importantly, the code is having difficulty at r=0 because of the 1/r term. Try rlow=0.01 or something.

[@baacademy1166]

@@jakeblanchard thank you sir, it run. I have another question. how can I solve and plot for multiple equations. In other words, how can I solve the above equation with different values of B and T. Thank you

[@jakeblanchard]

@@baacademy1166 I have two comments. First, I played around with your code and it appears that when you choose rlow>0, the result is quite sensitive to what you choose for rlow. You are going to have to do some work to make sure your results are correct. Matlab does have a method for addressing the singular term directly. I tried it and Matlab was having a hard time converging. This is further indication that you need to pay more attention to the solution around r=0.

[@jejeejajar]

@jakeblanchard Thank you so much for your response. I'll try to solve it. I doubt why I missed your video, it's wonderful ;)

[@LucasGuanes]

Thank You!

[@E1337Si]

Thank you good sir.

[@anujoybardhan]

Thanks a lot for the explanation. I am trying to solve a case where I have a first order Diff eqn and I have the boundary value known only for the end point. dp/dt= α(x-β) where α and β are constants. I know x at t0 and tf but I know p only for tf. Can I still use this bvp4c to solve my problem?

[@jakeblanchard]

Do you have a system of first order equations, or just one? If it's just one, then you can't impose all of your boundary conditions.

[@anujoybardhan]

I have a 2 first order equations and 2 variables p1 and p2. The DEs are dp1/dt=0 and dp2/dt=α(x-β). Now the first one is quite straightforward since p1 is constant and its value is known at t=tf, hence solved. Now for the second variable p2, the first order derivative dp2/dt is known at all points from t=t0 to t=tf but the value of p2 is known only at t=tf. I have solved this part manually, however , I need to solve this in matlab/simulink since my rest of simulation is on the same.

[@jakeblanchard]

@@anujoybardhan OK, so there is still an issue. With two first order equations you can only have two boundary conditions. You'll have to eliminate one. Once you do that, bvp4c should work.

[@anujoybardhan]

@@jakeblanchard Just to make it clear , I have 2 equations dp1/dt=0 and dp2/dt=α(x-β) and I have 2 boundary conditions i.e. p1(t=tf) and p2(t=tf). So with this information can I use bvp4c. Additionally in the step where bvp4ode(x,y) was defined, can I have 3 inputs ? For me those will be bvp4ode(t,x,p) ?

[@jakeblanchard]

@@anujoybardhan I'm still confused. What is x? Perhaps you could tell me something about the problem you are trying to solve.

[@sabeeqaali3428]

Sir can you please make a tutorial which guide us how to use bvp5c to solve coupled ordinary differential equations. Regards

[@jakeblanchard]

I've never used bvp5c. As far as I can tell, the usage is exactly like bvp4c.

[@sabeeqabatool7587]

OK, but want to know the difference between bvp4c and bvp5c. Can you please guide me about the method used behind bvp4c or limitations of bvp4c? So that I can compare it with bvp5c and find the solution of my problem, sir +Jake Blanchard ?

[@aagproyecto1794]

thank you for everything!!!! XD

[@upendrabajpai2448]

i am getting this error :( Error in bvp4 (line 5) sol=bvp4c(@bvp4ode,@bvp4bc,solinit); in program function bvp4 xlow=0; xhigh=80; solinit=bvpinit([0:10:80],[1 -1]); sol=bvp4c(@bvp4ode,@bvp4bc,solinit); xint=linespace(xlow,xhigh); Sxint=deval(sol,xint); plot(xint,Sxint(1,:)) end function dydx=bvp4ode(x,y) dydx=[y(2);-428.57*y(1)]; end function res=bvp4bc(ya,yb) res=[ya(1);-100 yb(2)]; end any idea how to correct it?

[@jakeblanchard]

You have the following statement in your code: res=[ya(1);-100 yb(2)]; You can't use 100 yb(2). You need an operator in there.

[@upendrabajpai2448]

sir ,thank you so much for replying , i corrected the mistake but still it is showing the same error. program function bvp4 xlow=0; xhigh=80; solinit=bvpinit([0:10:80],[1 -1]); sol=bvp4c(@bvp4ode,@bvp4bc,solinit); xint=linespace(xlow,xhigh); Sxint=deval(sol,xint); plot(xint,Sxint(1,:)) end function dydx=bvp4ode(x,y) dydx=[y(2);-428.57*y(1)]; end function res=bvp4bc(ya,yb) res=[ya(1)-100;yb(2)]; end

[@jakeblanchard]

You spelled linspace as linespace, but otherwise this looks ok. This was an error you should have been able to figure out on your own.

[@jejeejajar]

Great Video!!! Thank you so much. Your explanation is very clear. Could you give me a suggestion, please. If my boundary condition is in form y'... for example ... y'(0) = 0 and y'(10) = 0. How I can solve it in MATLAB? I try to guess initial point of y(0), but it seems wrong, Thank you in advance ;)

[@nbh9701]

Thank you for making this video.

[@harelchezi9239]

.Thank you

[@jakeblanchard]

I'm not sure what the problem is. If you send me your script, I'd be happy to identify the source of the error messages. My email is blanchard "at" engr "dot" wisc "dot" edu

[@safdar9688]

function newpaperssss Ha=1; Nt=0.5; s=2; Sc=1; Nb=0.5; Pe=1; Pr=1; Le=1; Nr=0.5; Nc=0.5; lambda=1; sigma=0.1; gamma=0.1; delta=0.5; solinit=bvpinit(linspace(0,10,100),[0 0 0 0 0 0 0 0 0]); sol=bvp4c(@bvpexam2,@bcexam2,solinit); x1=sol1.x; y1=sol1.y; %----------------------------------------------------- figure(1) plot(x1,y1(2,:),'b'); hold on figure(2) plot(x1,y1(4,:)); hold on figure(3) plot(x1,y1(6,:)); hold on figure(4) plot(x1,y1(8,:)); hold on z=y1(3,1); p=y1(7,1); function res = bcexam2(y0, yinf) res=[y0(1)-s;y0(2)-lambda;y0(4)-1;y0(6)-1;y0(8)-1;yinf(2)-0;yinf(4)-0;yinf(6)-0;yinf(8)-0]; end function ysol = bvpexam2(~,y) yy1 = (1/(1-D*y(1)^2))*(Ha*y(2)+y(2)^2-lambda*(y(4)-Nr*y(6)-Nc*y(8))-2*D*y(1)*y(2)*y(3)-(y(1)*y(3)+D*Ha*y(1)*y(3))); yy2 = (-3/(3+4*Rd))*(Pr*y(1)*y(5)+Pr*Nt*y(5)^2+Pr*Nb*y(5)*y(7)+Pr*delta*y(4)); yy3 = (Le*gamma*y(6)-Le*y(1)*y(7)-(Nt/Nb)*yy2); yy4 = (Pe*y(7)*y(8)+Pe*yy3*(sigma+y(9))-Sc*y(1)*y(9)); ysol = [y(2);y(3);yy1;y(5);yy2;y(7);yy3;y(9);yy4]; end end

[@jakeblanchard]

The two functions at the bottom should be after the main function. So the "end" at the bottom should be moved to the line after "p=y1(7,1);" Also, you should define the constants, like Ha in the function where they are used. There are ways to define them in the main routine and use them in the other functions, but for now I would suggest you define them where you use them. Also, you never defined D.

[@safdar9688]

@@jakeblanchard please give me WhatsApp number

[@safdar9688]

@@jakeblanchard i need help

[@safdar9688]

@@jakeblanchard sir I defined D and other parameters but cannot run please sir help kindly share WhatsApp number or Gmail, I send this file I cannot remove the error last 3 months please help me, sir, i need your help