How I can give condition & plot the solution of this differential equation. . . . . . . Please Guide

This is the equation for which

boundery condition are

theta(z=0)=0 degree

theta(z=h)=90 degree

where h=6

z=0:h

how to give condition here

e=8.85*10^-12
dele=11
E=1
k11=9
k33=9
k22=11
syms theta(z) z dtheta
dtheta=diff(theta,z)
d2theta=diff(theta,z,2)
eqn=d2theta+((k33-k11)*cos(theta)*sin(theta))*(dtheta)^2*(1/(k11*(cos(theta))^2+k22*(sin(theta))^2))+e*dele*E^2*cos(theta)*sin(theta)*(1/(k11*(cos(theta))^2+k22*(sin(theta))^2))
cond(theta(0)==0, theta(pi/2)==0)
thetaSol = dsolve(eqn,cond)
thetaSol = unique(simplify(thetaSol))
fplot(thetaSol)

[Torsten]

Neither equation nor boundary conditions of your mathematical description and your code coincide.

So the question: which one is correct ?

[Sam Chak]

@DEEPAK KARARWAL, So you want to find the value of A and initial value such that the simulation runs from to subject to , and ?

[Star Strider]

The ‘cond’ assignment was definitely a problem.

That corrected, the equation appears to be uniformly zero.

e=8.85*10^-12

e = 8.8500e-12

dele=11

dele = 11

E=1

E = 1

k11=9

k11 = 9

k33=9

k33 = 9

k22=11

k22 = 11

syms theta(z) z dtheta
dtheta=diff(theta,z)

dtheta(z) = 

d2theta=diff(theta,z,2)

d2theta(z) = 

eqn=d2theta+((k33-k11)*cos(theta)*sin(theta))*(dtheta)^2*(1/(k11*(cos(theta))^2+k22*(sin(theta))^2))+e*dele*E^2*cos(theta)*sin(theta)*(1/(k11*(cos(theta))^2+k22*(sin(theta))^2))

eqn(z) = 

cond = [theta(0)==0, theta(pi/2)==0];
thetaSol = dsolve(eqn,cond)

thetaSol = 

0

thetaSol = unique(simplify(thetaSol))

thetaSol = 

0

fplot(thetaSol, [0 pi/2])

.

[Torsten]

@Star Strider

I think the OP copied the code from somewhere in the internet and wants to use it to define his own problem as stated in the mathematical description.

[DEEPAK KARARWAL]

I want to plot theta as a function of z, if we have a liquid crystal cell, where cell thickness is h=6 and a gradient z varies from 0 to 6, so i want to plot z(on x axis), theta(on y axis),

Value of A=sqrt(dele*epsilnot*E^2/k11) for above written equation in main panel.

Original Equation is (k11cos^2(theta)+k33*sin^2(theta))d2theta+(k33-k11)*cos(theta)*sin(theta)*(dtheta)^2+(sqrt(e*epsilnot*E^2)*sin(theta)*cos(theta)

This is what i have simplified and make a code written above.

so finally boundry condition are:

theta(z=0)=0

theta(z=h)=pi/2

Parameters are:

epsilnot=8.85*10^-12

dele=11

E=1

k11=9

k33=9

k22=11

[Sam Chak]

What is the value of epsilnot? @DEEPAK KARARWAL

[DEEPAK KARARWAL]

value of epsilnot=8.85*10^-12

Answer 1

[Torsten]

dsolve doesn't succeed. Thus use a numerical solver (bvp4c) to solve your equation.

syms A theta(z) 
dtheta=diff(theta,z)

dtheta(z) = 

d2theta=diff(theta,z,2)

d2theta(z) = 

eqn = d2theta + A/2*sin(2*theta)==0;
cond = [theta(0)==0, theta(6)==pi/2];
thetaSol = dsolve(eqn,cond)

Warning: Unable to find symbolic solution.

 
thetaSol =
 
[ empty sym ]
 

[Sam Chak]

Can OP use the shooting method here (since the pendulum-like ODE is relatively simple) without the bvp4c?

[Torsten]

@Sam Chak

Do you think that the shooting method is easier than using bvp4c ?

epsilnot = 8.85*10^-12;
dele     = 11;
E        = 1;
k11      = 9;
k33      = 9;
k22      = 11;
A        = sqrt(dele*epsilnot*E^2/k11);
f = @(t, x) [x(2); ...
             - (A/2)*sin(2*x(1))];
tspan = [0 6];

xdot0 = 1.0;
xdot = fzero(@(xdot)fun(xdot,f,tspan),xdot0)

xdot = 0.2618

initc = [0 xdot];	% initial condition
[t, x] = ode45(f, tspan, initc);
plot(t,x(:,1))

function res = fun(xdot,f,tspan)
  initc = [0 xdot];	% initial condition
  [t, x] = ode45(f, tspan, initc);
  res = x(end,1) - pi/2;
end

[Sam Chak]

@Torsten

I understand bvp4c is a built-in BVP solver (collocation technique) that requires three ingredients to be specified: (1) the ODEs, (2) the boundary conditions and (3) the initial solution guess. My Professor only taught the shooting method because bvp4c was unavailable at the time of the Numerical Method course. In OP's case, the problem should be simple enough for him to learn about the shooting method.

[DEEPAK KARARWAL]

okay sir, if finally in place of pi/2, i want to replace asind((E-ec)/ec)

where E=anything ,ec=1.2 , but while i am giving this it is showing error

dele = 11 ;

E =5 ;

ec=1.2;

k11 = 9*10^-12;

k33 = 10*10^-12 ;

k22 = 11*10^-12 ;

epsilnot=8.85*10^-12;

A = sqrt(dele*epsilnot*E^2/k11 );

f = @(t, x) [x(2 ); ...

- (A/2)*sind(2*x(1 ))];

tspan = [0 6];

xdot0 = 1.0 ;

xdot = fzero(@(xdot)fun(xdot,f,tspan),xdot0);

initc = [0 xdot]; % initial condition

[t, x] = ode45(f, tspan, initc );

plot(t,x(:,1 ))

hold on

function res = fun(xdot,f,tspan )

initc = [0 xdot]; % initial condition

[t, x] = ode45(f, tspan, initc );

res = x(end,1) - asind((E-ec)/ec);

end

Error in abcd_final (line 14)

xdot = fzero(@(xdot)fun(xdot,f,tspan),xdot0);

where should I change this?

[Torsten]

For asind to be real-valued, the argument must be between -1 and 1. Since this is not the case for your settings, you get a complex number as result.

Do you want to set a complex number as boundary value at t=6 ? I guess no.

dele = 11 ;
E =5 ;
ec=1.2;
bv = asind((E-ec)/ec)

bv = 9.0000e+01 - 1.0427e+02i

k11 = 9*10^-12;
k33 = 10*10^-12 ;
k22 = 11*10^-12 ;
epsilnot=8.85*10^-12;
A = sqrt(dele*epsilnot*E^2/k11 );
f = @(t, x) [x(2 ); ...
    - (A/2)*sind(2*x(1 ))];
tspan = [0 6];

xdot0 = 1.0 ;
xdot = fzero(@(xdot)fun(xdot,f,tspan,bv),xdot0);

Error using fzero
Initial function value must be finite and real.

initc = [0 xdot]; % initial condition
[t, x] = ode45(f, tspan, initc );
plot(t,x(:,1 ))
hold on

function res = fun(xdot,f,tspan,bv )
initc = [0 xdot]; % initial condition
[t, x] = ode45(f, tspan, initc );
res = x(end,1) - bv;
end

[DEEPAK KARARWAL]

yes sir, but what I want, is to give such a boundery condition at theta(z)=(__) such that I will get varying angle from 0 to 90 degree as we increase the value of E, where E is electric field and contained in the expression of A.

[Torsten]

Not clear what you mean.

The boundary value for theta at z = 6 can be set by writing it in the variable "bv" of my code above. Experiment with it.

Answer 2

[Sam Chak]

Hi @DEEPAK KARARWAL

Giiven the parameters, it seems that if you select initial values and , the boundary values are satisfied.

epsilnot = 8.85*10^-12;
dele     = 11;
E        = 1;
k11      = 9;
k33      = 9;
k22      = 11;
A        = sqrt(dele*epsilnot*E^2/k11);

f = @(t, x) [x(2); ...
	        - (A/2)*sin(2*x(1))];

tspan = [0 6];

initc = [0 pi/12];	% initial condition

[t, x] = ode45(f, tspan, initc);
plot(t, x(:,1), 'linewidth', 1.5), grid on, xlabel('t'), ylabel('\theta')

x(end,1)    % π/2 at θ(6)

ans = 1.5708

[DEEPAK KARARWAL]

Thanks you sir , I will analyze it and than report you again, Thank you

Answer 3

[MOSLI KARIM]

 %%

function answer 
clc
clear all
close all
global A 
epsilnot = 8.85*10^-12;
dele     = 11;
E        = 1;
k11      = 9;
k33      = 9;
k22      = 11;
A        = sqrt(dele*epsilnot*E^2/k11);

solinit=bvpinit(linspace(0,6),[0;pi/12])
sol=bvp4c(@fct,@bc,solinit)
figure(1)
plot(sol.x,sol.y(1,:))

    function  dxdy=fct(x,y)
    dxdy=[y(2); -(A/2)*sin(2*y(1))];   
    end 

    function res=bc(ya,yb)
        
    res=[ya(1);yb(1)-90]    
    end 

end